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Three-dimensional  curvilinear-grid  hydrodynamic  models  for  partially-mixed  es- 
tuaries developed  by  Sheng  (1989)  and  Sheng  et  al.  (1989)  were  significantly  en- 
hanced, applied  to  the  James  River,  and  tested  to  determine  the  sensitivity  of  model 
results  to  numerous  model  parameters.  The  model  computes  the  inter-tidal  and 
longterm  dynamics  of  surface  elevation  and  three-dimensional  field  of  velocity,  and 
salinity  in  partially-mixed  estuaries  with  complex  geometry  and  bathymetry. 

The  momentum  and  continuity  equations  in  the  curvilinear  coordinates  are  solved 
on  a space-staggered  grid  system  using  an  implicit  finite-difference  algorithm  for  the 
“external”  vertically-integrated  flow  and  vertically  implicit  procedure  for  the  “inter- 
nal” flow.  Details  of  the  numerical  procedures  are  presented.  The  numerical  accuracy 
of  the  model  is  verified  by  comparing  model  results  with  analytic  solutions  for  a vari- 
ety of  idealized  problems  describing  circulations  in  simple  basins  forced  by  constant 
surface  slope,  seiche,  tide,  wind,  and  density  gradient. 

Model  application  to  the  James  River  reproduces  the  observed  density-driven 
current.  Effect  of  spring-to-neap  tidal  dynamics  on  vertical  stratification  and  density- 
driven  currents  is  simulated  by  the  model.  Longitudinal  position  of  the  head  of 
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salinity  intrusion  was  found  to  move  up  and  down  the  estuary  due  to  change  in  river 
discharge  and  tidal  condition.  Increasing  river  discharge  caused  an  increase  in  relative 
stratification.  Spring  tide  produced  more  mixing  and  enhanced  salt  intrusion  while 
neap  tide  produced  more  stratification.  The  observed  frontal  structure  near  Newport 
News  Point  was  successfully  simulated  as  a narrow  zone  of  surface  flow  convergence 
and  strong  vertical  transport.  Simulation  of  the  formation  and  evaluation  of  the  front 
show  reasonable  agreement  with  available  field  data.  Sensitivity  of  the  model  results 
to  such  parameters  as  turbulence  model,  advection  scheme,  horizontal  grid,  vertical 
grid,  and  bottom  boundary  condition  was  investigated  and  the  results  were  presented. 


xviii 


CHAPTER  1 
INTRODUCTION 


1.1  Background 

Estuaries  are  the  transition  regions  where  the  fresh  water  from  the  land  mixes 
with  the  saline  water  from  the  sea.  Cameron  and  Pritchard  (1963)  have  defined 
an  estuary  as  “a  semi-enclosed  coastal  body  of  water  which  has  a free  connection 
with  the  open  sea  and  within  which  sea  water  is  measurably  diluted  with  fresh  water 
derived  from  land  drainage”  (p.  306).  This  definition  implies  two  important  dynamics 
in  an  estuary:  a horizontal  density  gradient  between  the  fresh  water  and  the  sea 
water,  and  a generation  of  kinetic  energy  caused  by  tidal  currents.  These  dynamics 
control  the  circulation  pattern,  the  mixing  processes,  and  the  salinity  and  temperature 
distributions  that  are  the  major  physical  problems  that  need  to  be  understood  in  an 
estuary.  Ecologically,  most  estuaries  are  controlled  by  physical  factors  (Kinne,  1967) 
such  as  tide,  river  discharge,  wind,  salinity,  and  morphology  of  the  estuarine  area. 

Vertical  mixing  caused  by  the  bottom  and/or  surface  friction,  tidal  currents,  and 
the  river  discharge  produces  vertical  and  longitudinal  gradients  of  salinity  in  the 
mixing  zone.  Large  river  discharge  and  small  mixing  energy  form  sharp  vertical 
gradients  in  an  estuary,  which  is  referred  to  as  “stratified”.  Salt-wedge  is  the  limiting 
case  of  the  stratified  estuary.  If  the  river  discharge  is  small  and  the  tidal  mixing 
energy  large,  slight  vertical  gradients  of  salinity  will  be  developed.  In  the  limiting 
case,  a “well-mixed”  estuary  is  formed. 

Partially-mixed  estuaries  are  intermediate  between  salt-wedge  and  well-mixed  es- 
tuaries. James  River  is  a partially- mixed,  coastal  plain  type  estuary  like  the  majority 
of  estuaries  along  the  east  coast  of  the  United  States.  These  estuaries  are  charac- 
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terized  dynamically  by  vertical  salinity  variations  of  1-10  °/00  (Bowden,  1967)  and 
a net  circulation  pattern  which  moves  seaward  in  the  upper  layers  and  landward  in 
the  lower  layers  due  to  the  pressure  gradient  caused  by  the  density  difference  of  the 
salt  and  fresh  waters.  Rattray  and  Hansen  (1962)  showed  that  residual  currents  were 
directly  proportional  to  the  longitudinal  density  gradient  and  inversely  proportional 
to  the  vertical  eddy  viscosity.  Stratification  exists  in  partially-mixed  estuaries,  and 
this  changes  in  various  ways  according  to  the  relative  amount  of  river  discharge  and 
available  tidal  mixing  energy.  If  the  mixing  energy  is  cyclical,  like  a spring-neap  tidal 
variation,  alternate  stratification  and  destratification  may  be  evident  (Haas,  1977; 
Cerco,  1982).  Current-induced  turbulent  mixing  plays  an  important  role  in  affect- 
ing the  stratification  and  destratification  in  a partially-mixed  estuary.  Vertical  salt 
flux  occurs  through  vertical  advection  and  entrainment  of  salinity  across  the  interface 
between  fresh  water  and  salt  water.  Tidal  motion  enhances  the  vertical  turbulent  ex- 
change of  momentum  and  salt.  A two-layered  flow  exists  in  a partially-mixed  estuary 
with  the  surface  of  no  net  motion  near  the  maximum  salinity  gradient  (halocline).  In 
addition  to  the  longitudinal  and  vertical  salinity  gradients,  a lateral  salinity  gradient 
also  occurs  because  of  the  Coriolis  force  and  complex  shorelines.  Thus,  motions  in 
partially-mixed  estuaries  are  basically  three-dimensional  and  time-dependent  because 
of  irregular  shape  of  the  basin,  density-induced  baroclinicity  and  the  Coriolis  force. 

Salinity  and  temperature  influence  the  density  of  estuarine  waters,  but  salinity  is 
more  effective  than  temperature  because  the  range  of  salinity  variation  is  much  larger 
than  that  of  temperature  variation.  Salinity  ranges  from  freshwater  to  saltwater  in 
the  estuary.  Wind  also  influences  the  vertical  salinity  and  velocity  structures  while 
interacting  with  internal  friction  and  bottom  friction.  Surface  currents  are  generally 
in  the  same  direction  as  the  local  wind  stress  even  though  they  are  affected  by  the 
Coriolis  force.  Thus  downestuary  winds  increase  stratification  while  upestuary  winds 
decrease  it.  Nonlocal  wind  effects  can  affect  stratification  through  change  of  sea  level 
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at  the  mouth  of  the  estuary  (Wang  and  Elliott,  1978). 

River  or  open  sea  boundary  conditions  vary  with  time,  and  the  location  of  salinity 
intrusion  changes  accordingly.  If  the  head  of  salt  intrusion  is  defined  as  the  point 
where  freshwater  meets  saltwater,  the  head  of  salt  intrusion  is  dependent  on  the 
amount  of  river  discharge,  the  salinity  at  the  open  sea  boundary,  tidal  currents,  wind 
stresses,  and  the  geometry  of  the  basin.  They  are  all  time-dependent  except  the 
geometry  of  the  basin.  Accurate  specification  of  the  boundary  conditions  is  essential 
to  accurate  prediction  of  circulation  and  salinity  distribution. 

Turbulence  models  have  been  developed  to  simulate  the  mean  and  turbulent  mo- 
tions in  most  stratified  flow  situations.  A turbulence  model  is  defined  as  a system 
of  equations  which  calculate  Reynolds  stress  terms  in  the  mean-flow  equations  and 
similar  terms  in  temperature  and  salinity  equations.  They  parameterize  the  effect  of 
turbulence  on  the  mean  flow.  Through  the  Reynolds  decomposition  of  the  governing 
equations  and  Reynolds  averaging,  Reynolds  stress  equations  are  developed.  However, 
the  system  of  equations  is  not  closed  since  every  equation  contains  terms  consisting 
of  higher  order  correlations  which  are  unknown.  The  system  can  only  be  closed  by 
using  a semi-empirical  (e.g.  various  levels  of  second-order  closure  model)  or  ad-hoc 
empirical  models  (e.g.  simple  eddy-viscosity  models).  An  equilibrium  closure  model 
(Sheng  et  al.,  1989)  is  included  in  the  three-dimensional  model  used  for  this  study. 
To  reduce  computational  effort,  the  system  of  equations  can  be  simplified  with  such 
assumptions  as  local  equilibrium  of  turbulence.  Although  Reynolds  stress  models  can 
provide  rather  accurate  results  (e.g.  Sheng,  1984;  Sheng  and  Villaret,  1989),  Reynolds 
stress  models  are  quite  complicated,  and  require  substantial  computer  resources. 

At  present,  hydrodynamic  models  of  circulation  are  important  tools  used  by  both 
engineers  and  scientists  attempting  to  quantify  and  understand  complex  physical, 
chemical  and  biological  processes  in  estuaries  and  coastal  areas.  Historically,  engineers 
and  scientists  preferred  physical  modeling  to  mathematical  modeling  for  studying  cir- 
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culation  and  transport  in  most  of  the  complicated  environmental  systems  because 
of  difficulties  in  finding  mathematical  solutions.  However,  physical  modeling  has  in- 
herent scaling  problems  and  is  very  costly.  As  computers  become  faster  and  more 
accurate,  mathematical  models  have  replaced  physical  modeling  as  the  major  tool  for 
studying  circulation  and  transport.  Some  differences  exist  between  numerical  results 
and  analytical  solutions  of  mathematical  problems  because  of  the  digital  nature  of 
the  computations.  To  reduce  these  differences  and  enhance  efficiency  of  the  compu- 
tations, numerous  numerical  schemes  and  solution  techniques  have  been  developed. 
The  numerical  models  developed  must  be  checked  against  analytical  solutions  before 
applying  them.  Agreement  between  numerical  results  and  analytical  solutions  sug- 
gests accurate  and  consistent  numerics  of  the  numerical  model.  To  test  the  model’s 
ability  to  simulate  such  physical  processes  as  turbulent  mixing  and  bottom  friction,  it 
is  necessary  to  calibrate  and  verify  the  numerical  model  with  laboratory  and/or  field 
data  on  circulation  and  transport. 

1.2  Review  of  Previous  Works 

Modeling  of  circulation  and  salinity  transport  dynamics  in  the  James  River  has 
been  active  over  the  last  decade.  Circulation  modeling  has  included  the  flow  induced 
by  tide,  wind,  and  density  gradients,  with  the  primary  emphasis  on  the  tidally  driven 
circulation.  Cerco  (1982)  applied  a two-dimensional  laterally  averaged  model  to  the 
James  River  to  predict  intertidal  variation  of  surface  level,  velocity,  and  salinity. 
This  model  could  simulate  the  destratification-stratification  cycle  coincident  with  the 
spring-neap  tidal  cycle  in  the  James  River  Estuary.  Heltzel  and  Granat  (1988)  applied 
a two-dimensional  vertically-averaged  model  to  the  lower  James  River  to  study  the 
impact  of  artificial  islands  for  dredge  spoil  disposal,  channel  deepening,  and  port 
expansion.  This  model  could  simulate  the  convergence  phenomena  in  the  tidal  flow, 
but  could  not  show  the  complex  vertical  structure  of  the  flow  field. 

Most  of  the  early  work  on  the  numerical  modeling  of  estuaries  used  vertically- 
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integrated  equations  (Hansen,  1956;  Leendertse,  1967;  Blumberg,  1975;  Wang,  1975; 
Johnson,  1980).  This  approach  was  found  to  be  useful  for  studying  tidal  dynamics 
and  storm  surges  although  it  cannot  simulate  the  internal  features  of  the  estuar- 
ine circulation.  In  other  cases,  laterally  averaged  equations  were  used  to  study  the 
vertical  structure  of  currents  and  species  (Spaulding,  1974;  Blumberg,  1976;  Cerco, 
1982),  although  no  lateral  flow  structure  could  be  provided.  The  development  of 
three-dimensional  numerical  models  started  in  the  late  1960’s,  but  have  progressed 
significantly  during  the  past  decade  because  of  the  availability  of  high  speed  comput- 
ers. Details  of  a few  three-dimensional  numerical  estuarine  models  can  be  found  in 
Caponi  (1976),  Sheng  (1983,1987),  Lynch  (1983),  Blumberg  and  Mellor  (1983),  and 
Oey  et  al.  (1985). 

The  early  three-dimensional  models  were  developed  by  dividing  the  vertical  water- 
column  into  layers  and  solving  each  layer  with  a vertically  averaged  two-dimensional 
model  (e.g.,  Laevastu,  1974).  Laevastu  (1974)  introduced  a two-layer  model  consist- 
ing of  a stack  of  two  vertically  averaged  models  and  a staggered  finite  difference  grid 
system.  Wang  and  Connor  (1975)  also  used  the  layer  concept  using  a finite  element 
method.  This  layered  approach  has  the  disadvantage  that  a layer  thickness  can  disap- 
pear in  simulations  with  strong  vertical  velocities  such  as  downwelling  and  upwelling. 
To  alleviate  this  weakness,  many  researchers  have  used  the  leveled  approach.  In  a 
leveled  model,  the  level  interfaces  are  fixed,  and  free  exchanges  of  water  mass,  mo- 
mentum, salt  and  heat  are  allowed  across  the  interfaces.  The  leveled  approach  is  used 
by  Leendertse  and  Liu  (1975),  Caponi  (1976),  and  Sheng  (1983,  1987). 

Sheng  (1986c,  1986d,  1987)  developed  a three-dimensional  hydrodynamic  model 
which  can  be  used  in  non-orthogonal  boundary-fitted  grid  as  well  as  orthogonal  or 
Cartesian  grid.  As  in  the  earlier  Cartesian  grid  model  (Sheng,  1983),  the  curvilinear- 
grid  model  uses  similar  model-splitting  technique,  vertical  turbulence  closure  scheme, 
modified  alternating  direction  implicit  (ADI)  scheme  for  the  external  mode,  and 
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vertically-implicit  numerical  schemes.  The  governing  equations  in  the  curvilinear  grid 
system  were  derived  by  tensor  transformation  in  terms  of  the  contravariant  velocity 
components. 

Swanson  (1987)  developed  a three-dimensional  boundary-fitted  hydrodynamic 
model  in  terms  of  Cartesian  velocity  components.  He  developed  a Helmholtz  equation 
for  surface  elevations,  which  is  derived  by  eliminating  the  velocities  from  the  continu- 
ity equation  and  the  momentum  equations,  and  solved  it  with  the  Gaussian  elimina- 
tion or  successive  over  relaxation  (SOR)  method.  This  method  has  an  advantage  in 
handling  boundary  conditions  since  velocity  components  are  known  simultaneously 
at  the  new  time  step  compared  with  an  ADI  method.  However,  the  advantage  of  an 
ADI  scheme  is  far  greater  in  terms  of  its  efficiency.  The  tridiagonal  matrix  associated 
with  the  ADI  method  is  very  efficient  to  solve  on  any  computer,  while  SOR  needs 
more  computing  time  because  of  the  repeated  iterations  required.  Swanson  used  con- 
stant vertical  eddy  coefficients,  which  is  incorrect  (turbulence  is  time-dependent  and 
spatially-varying)  and  is  difficult  to  adjust  in  each  new  model. 

The  development  of  turbulence  closure  models  has  been  of  major  significance  be- 
cause of  the  need  to  accurately  simulate  the  vertical  fluxes  of  salt  and  momentum 
and  their  dependence  on  stratification.  A number  of  turbulence  models  have  been  de- 
veloped. For  example,  Donaldson  (1973)  and  his  colleagues  at  Aeronautical  Research 
Associates  of  Princeton  developed  a hierachy  of  turbulence  models.  Sheng  (1982, 
1984)  modified  the  Reynolds  stress  model  and  Sheng  and  Chiu  (1986)  and  Sheng 
and  Villaret  (1989)  modified  the  equilibrium  closure  and  turbulent  kinetic  energy 
(TKE)  closure  models  for  marine  and  freshwater  applications.  Mellor  and  Yamada 
(1982)  described  a similar  turbulence  model  which  was  applied  by  Blumberg  and 
Mellor  (1983)  and  Oey  et  al.  (1985).  They  classified  the  turbulence  model  as  level 
4 (Reynolds  stress),  level  3 (TKE  closure),  level  2 (equilibrium  closure),  and  level 
2.5  models  according  to  the  complexity  of  the  turbulence  equations.  The  turbulence 
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closure  models  have  been  applied  to  simulate  complex  mixing  processes  (Sheng,  1982, 
1985;  Blumberg  and  Mellor,  1983;  Oey  et  al,  1985).  In  addition  to  the  Reynolds  stress 
equations  of  turbulent  motion,  equations  for  turbulent  kinetic  energy  and  turbulent 
macroscale  had  to  be  added  to  the  system  of  equations  to  complete  the  turbulence 
model.  Although  various  levels  of  turbulence  closure  models  have  been  used  by  vari- 
ous modelers,  little  has  been  done  to  compare  the  results  obtained  with  various  levels 
of  closure  models  versus  realistic  data  from  estuaries. 

Conventional  finite-difference  models  use  only  Cartesian-grid  and  have  shortcom- 
ings for  applications  where  physical  boundaries  are  very  complex  and  thus  may  require 
excessively  fine  grid  in  order  to  resolve  the  boundary  in  the  model.  To  allow  a more 
flexible  arrangement  and  spatial  variation  of  the  grids,  one  can  use  curvilinear-grid 
finite-difference  models  which  were  discussed  above  or  finite-element  models  (Wang, 
1975;  Wang  and  White,  1976).  The  finite-element  models  allow  flexibility  in  the  dis- 
cretization of  the  spatial  domain.  However,  because  the  matrix  of  the  discretized 
system  is  irregular,  direct  inversion  method  which  requires  very  extensive  compu- 
tational effort,  is  usually  used.  In  the  finite-difference  models,  however,  tridiagonal 
matrix  system  generally  exists  and  can  be  efficiently  solved.  The  curvilinear-grid 
finite-difference  models  (e.g.  Spaulding,  1984;  Sheng,  1986c,  1986d,  1987).  This  al- 
low extensive  user  flexibility  in  spacing  of  the  horizontal  grid  structure  while  still 
maintaining  a rather  simple  matrix  system  and  allowing  efficient  solutions.  However, 
there  are  some  disadvantages  of  this  approach  such  as  more  computational  efforts 
due  to  the  extra  terms  in  the  transformed  equations  and  potential  numerical  dissipa- 
tion/dispersion problems  caused  by  skewness  and  rapid  change  in  grids. 

1.3  Modeling  Strategy 

The  examination  and  analysis  of  available  data,  such  as  river  flow,  water  eleva- 
tions, wind,  current,  and  salinity,  should  be  performed  to  classify  the  dynamics  of 
the  water  body  along  with  its  range  of  variability.  Models  that  do  not  adequately 
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resolve  the  depths  of  navigation  channels  of  an  estuarine  system  may  have  difficulty 
in  predicting  the  correct  degree  of  stratification  and  the  strength  of  tidally  averaged, 
near-bottom  upstream  flow.  The  Hansen  and  Rattray  theory  (1965)  indicates  the 
importance  of  tidal  dispersion  in  maintaining  the  salt  balance  in  an  estuary. 

The  tide  in  an  estuary  may  have  the  characteristics  of  a progressive,  standing,  or 
mixed  wave.  Flood  and  ebb  tide  are  not  only  unequal  in  velocity  at  a fixed  point  in 
most  estuaries,  but  also  do  not  always  flow  opposite  to  each  other.  Ludwick  (1975) 
attributed  this  inequality  in  tidal  flow  to  three  factors:  (1)  mixing  between  river 
water  and  sea  water  produced  by  turbulent  exchange  and  vertical  entrainment  across 
a density  interface,  (2)  the  Coriolis  force  which  deflects  the  ebb  and  flood  dominant 
flows  along  different  perimeters  of  the  estuary,  and  (3)  evasion  of  ebb-  and  flood-tidal 
pathway.  Although  all  three  factors  may  occur  concurrently,  the  effect  of  the  Coriolis 
force  is  usually  overshadowed  by  the  other  two  factors. 

Exchange  of  various  forms  of  energy  produces  the  mixing  in  a vertically  stratified 
system.  For  example,  turbulent  kinetic  energy  extracted  from  the  mean  flow  in  a 
frictional  boundary  layer  produces  mixing  when  exchanged  with  potential  energy. 
Various  parameters  have  been  introduced  to  classify  the  degree  of  stratification.  Ippen 
(1966)  introduced  the  stratification  parameter  defined  as  the  ratio  of  rate  of  turbulent 
energy  dissipation  divided  by  rate  of  gain  of  potential  energy.  Hansen  and  Rattray 
(1965)  developed  the  dimensionless  parameter  by  which  the  degree  of  stratification 
can  be  related  to  buoyancy  and  mixing  energy. 

Festa  and  Hansen  (1976)  discussed  the  effects  of  depth  and  river  discharge  on 
estuarine  circulation  and  salinity  in  their  two-dimensional  model  based  on  vorticity 
and  salt-balance  equations.  They  found  that  decreasing  the  discharge  allowed  the 
head  of  salt  intrusion  to  move  upstream.  As  the  discharge  is  decreased,  estuarine 
circulation  is  also  weakened,  although  it  occurs  over  the  larger  area.  Greater  depths 
enhanced  circulation  and  moved  the  head  of  salt  intrusion  upstream. 
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Geyer  and  Farmer  (1989)  measured  the  time  variation  of  salt  intrusion  in  an  estu- 
ary during  a variety  of  tidal  and  run-off  conditions  using  high-resolution  instruments 
and  sampling  methods.  The  salt  wedge  was  found  to  vary  considerably  in  position 
and  vertical  structure  throughout  the  tidal  cycle  due  to  interaction  of  tidal  flow  and 
density-induced  motion.  During  the  flood  cycle,  salt  wedge  progresses  up  the  estuary 
as  a gravity  current,  while  during  the  ebb  the  salinity  structure  was  eroded  by  shear 
instability,  in  which  vertical  shears  become  so  large  that  the  pycnocline  becomes  un- 
stable and  causes  the  salt  wedge  to  break  down.  Here  the  most  important  dynamic 
variable  is  the  gradient  Richardson  number  (Ri),  which  reflects  the  stability  of  the 
shear  layer  and  thus  determines  whether  momentum  and  salt  can  be  exchanged  across 
the  pycnocline. 

Tides  in  estuaries  are  modulated  by  river-forced  estuarine  plumes.  A numerical 
study  of  this  phenomena  was  done  by  Chao  (1990).  In  a typical  midlatitude  estuary- 
shelf  environment,  the  river-forced  circulation  can  be  divided  into  three  regions:  a 
seaward  surface  flow  and  a landward  undercurrent  inside  the  estuary,  a strong  anti- 
cyclonic  surface  gyre  overlying  a weak  cyclon  off  the  estuary  mouth,  and  a coastal  jet 
flowing  in  the  direction  of  Kelvin  wave  propagation  farther  down-stream.  The  tidal 
residual  eddies  enhance  the  plume  growth  but  tend  to  trap  the  plume. 

Frontal  systems  in  estuaries  are  also  very  important  in  estuarine  circulation  and 
transport.  Limited  observations  by  Byrne  et  al.  (1987)  suggested  that  oyster  larvae 
remain  in  the  water  column  for  a period  of  two  to  three  weeks  and  are  carried  by  water 
movement  associated  with  tidal  front  and  tidal  residual  circulation.  A mathematical 
model  that  describes  the  formation  and  dilution  of  a frontally  bounded  river  plume  is 
presented  by  O’Donnell  (1990).  The  behavior  of  plumes  in  time-dependent  flows  like 
tidal  current  were  examined  in  his  numerical  test.  The  interaction  of  the  plume  front 
and  the  ambient  flow  was  found  to  have  major  influence  on  the  spreading  and  mixing 
of  the  buoyant  water.  River  plumes  may  only  mix  significantly  with  the  underlying 
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ambient  fluid  at  particular  phase  of  the  tide,  but  the  mixing  generally  occurs  over 
a large  area.  Thus  the  boundary  conditions  applied  in  estuarine  circulation  models 
may  need  to  be  carefully  reconsidered. 

In  any  case,  the  density  distribution  determines  the  stratification.  Therefore,  a 
synoptic  density-field  survey  should  be  made.  To  remove  tidal  effects  as  much  as 
possible,  it  is  desirable  to  sample  at  slack.  During  slack  water,  stratification  may  be 
stronger  than  during  the  other  stages,  making  water  column  sampling  easier  in  low 
currents.  The  usual  way  to  sample  at  slack  is  that  the  boat  is  following  slack  water. 

1.4  Scope  of  This  Study 

Modeling  of  a partially  mixed  estuary  is  one  of  the  most  challenging  areas  in  com- 
putational coastal  hydrodynamics.  Partial  stratification  takes  place  in  many  coastal 
plain  estuaries  like  the  James  River  (Figure  1.1),  where  freshwater  inflows  and  tidal 
mixing  interact  significantly.  Partial  stratification  requires  that  vertical  variation  of 
flow  variables  be  adequately  resolved  such  that  the  model  can  simulate  the  density- 
driven  currents.  Thus,  in  contrast  to  traditional  barotropic  tidal  hydrodynamics,  the 
simulation  of  baroclinic  circulation  requires  the  accurate  representation  of  vertical  ve- 
locity and  density  structures,  and  the  use  of  three-dimensional  models  (Sheng,  1983; 
Oey  et  al.).  Sheng  et  al  (1989a)  simulated  the  tidal,  wind-driven,  and  baroclinic  cir- 
culation in  Chesapeake  Bay,  while  Sheng  et  al  (1989b)  simulated  the  tidal  circulation 
and  baroclinic  circulation  in  the  James  River  using  a very  coarse  grid.  This  study 
provide  more  detailed  simulations  of  tidal  and  baroclinic  circulation  in  the  James 
River. 

Traditionally,  Cartesian  grids  have  been  used  for  the  estuarine  hydrodynamic 
model.  For  most  estuaries,  the  geometry  and  bathymetry  are  quite  complex  and  ir- 
regular so  that  the  generation  of  an  orthogonal  grid  is  very  difficult.  Models  developed 
by  Blumberg  and  Mellor  (1987)  require  an  orthogonal  grid  and  cannot  accept  a non- 
orthogonal  grid,  and  hence  cannot  accurately  resolve  the  complex  shorelines.  On  the 
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other  hand,  boundary-fitted  grid  hydrodynamic  models  appear  to  be  an  excellent  tool 
to  allow  the  user  to  accurately  resolve  complex  geometries  for  estuarine  and  coastal 
applications.  Few  studies  focused  on  model  development  and  validation  against  an- 
alytic solutions.  For  example,  Sheng  and  Hirsh  (1984)  developed  a two-dimensional 
curvilinear-grid  model  for  coastal  ocean  circulation.  Spaulding  (1984)  applied  a two- 
dimensional  model  to  simulate  tides  in  the  North  Sea.  Practical  applications  generally 
require  a fully  three-dimensional  model  analysis  because  of  stratification,  complex  ge- 
ometry, and  topography.  Although  the  three-dimensional  boundary-fitted  grid  hydro- 
dynamic  model  developed  by  Sheng  has  been  tested  against  some  analytic  solutions 
(Sheng  and  Hirsh,  1984;  Sheng,  1987),  this  study  provides  more  detailed  analytical 
tests  of  the  three-dimensional  model. 

Despite  the  complexity  of  the  Reynolds  stress  models,  such  models  do  contain 
more  physics  and  hence  the  model  constants  remain  unchanged  for  all  model  ap- 
plications (Sheng,  1982).  A Reynolds  stress  model  was  used  to  produce  successful 
simulations  of  the  wave  boundary  layer  (Sheng,  1982;  Sheng,  1984;  Sheng,  1986a; 
Sheng  and  Villaret,  1989),  wave-current  boundary  layer  (Sheng,  1984),  and  vegeta- 
tion boundary  layer  (Sheng,  1982).  However,  three-dimensional  models  generally  do 
not  use  Reynolds  stress  models  to  represent  the  turbulent  mixing,  because  of  the  com- 
plexity. Instead,  simplified  second-order  closure  models  have  been  developed  and  used 
with  three-dimensional  models.  These  include  the  turbulent  kinetic  energy  closure 
model  (Sheng  and  Villaret,  1989)  and  the  equilibrium  closure  model  (Sheng  and  Chiu, 
1986;  Sheng  et  al. , 1989a).  Both  of  these  models  contain  the  assumption  that  the 
turbulence  time  scale  is  much  less  than  the  mean  flow  time  scale.  The  equilibrium 
closure  model  contains  the  additional  assumption  that  turbulence  does  not  change 
significantly  over  the  turbulent  macroscale.  The  TKE  closure  model  resolves  the 
turbulent  dynamics  by  including  the  differential  transport  equation  for  q 2,  while  the 
other  second-order  correlation  equations  are  simplified  to  algebraic  equations.  The 
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TKE  closure  model  is  similar  to  the  Level-3  model  of  Mellor  and  Yamada  (1982)  and 
the  k — e model  of  Rodi  (1980),  but  uses  an  equation  for  the  macroscale  (A)  instead  of 
the  equation  for  turbulence  dissipation  (e)  used  by  the  other  two  models.  The  equilib- 
rium closure  model  is  slightly  more  simplified  than  the  TKE  closure  model,  since  only 
algebraic  equations  are  used  to  resolve  the  second-order  correlations.  The  equilibrium 
closure  model  is  similar  to  Mellor  and  Yamada’s  Level-2  model.  Instead  of  solving 
the  differential  transport  equations  for  A which  consists  completely  of  modeled  terms, 
simplified  second-order  closure  models  often  solve  a set  of  integral  constraints  for  A 
(Sheng  and  Chiu,  1986). 

Results  of  the  simplified  second-order  closure  models  have  been  quite  successful. 
Sheng  and  Villaret  (1989)  used  the  TKE  closure  model  to  simulate  a wave  boundary 
layer,  a thermally  stratified  boundary  layer  and  a sediment-laden  boundary  layer. 
On  the  other  hand,  the  equilibrium  closure  model  has  also  been  found  to  be  very 
successful  for  simulating  the  vertical  mixing  in  estuaries  (Sheng  et  al.,  1989b,  1989c; 
Johnson  et  al .,  1991),  despite  the  further  simplification.  There  exists  no  unequivocal 
evidence  to  suggest  that  the  TKE  closure  model  is  superior  to  the  equilibrium  closure 
model.  This  subject  will  be  investigated  in  this  study  by  comparing  results  of  TKE 
closure  model  and  equilibrium  closure  model  versus  field  data  from  the  James  River. 

Recent  field  studies  in  the  James  River  have  obtained  quantitative  or  qualitative 
information  on  the  temporal  and  spatial  variations  of  salinity  stratification,  frontal 
structure,  and  gyre  formation  (Bryne,  1987;  Hepworth  and  Kuo,  1989).  These  field 
data  provide  a good  opportunity  for  detailed  calibration  and  validation  of  the  three- 
dimensional  curvilinear  grid  hydrodynamic  model  developed  by  Sheng  (1989).  Model 
calibration  and  validation  will  be  conducted  in  terms  of  detailed  intertidal  dynamics, 
water  level  vertical  flow  and  salinity  structures,  spring-to-neap  and  residual  long-term 
circulation  and  salinity  transport,  and  formation  and  evolution  of  tidal  fronts  in  the 
Hampton  Roads  area.  Sheng  et  al.  (1989c)  reported  that  the  use  of  cr-grid  may  lead 
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to  numerically-induced  salinity  transport  from  the  navigation  channel  to  the  shallow 
shoal,  unless  special  care  is  taken  to  improve  the  evaluation  of  baroclinic  terms  and 
the  advection  terms.  This  subject  will  also  be  further  investigated  in  this  study. 

1.5  Objective  of  This  Study 

The  objectives  of  this  study  are  (1)  to  further  enhance  the  three-dimensional 
hydrodynamic  curvilinear- grid  model  developed  by  Sheng  (1989)  and  Sheng  et  al. 
(1989)  for  simulating  circulation  and  transport  patterns  in  coastal  areas  and  estu- 
arine systems,  (2)  to  apply  the  model  to  the  James  River  to  simulate  the  observed 
short-term  and  long-term  circulation  features,  (3)  to  determine  the  model  sensitiv- 
ity to  various  model  features  including  horizontal  grid  resolution,  bottom  friction, 
turbulence  closure  scheme,  and  vertical  grid. 

1.6  James  River 

The  Chesapeake  Bay  estuarine  system  was  formed  by  the  drowning  of  a cascading 
series  of  river  valleys,  the  most  seaward  of  which  was  the  James.  Estuaries  formed 
in  this  way  are  called  coastal  plain  estuaries.  In  this  type  of  estuary,  the  primary 
forcing  functions  for  salinity  and  current  distribution  are  tidal  mixing  which  depends 
on  the  tidal  range  and  basin  configuration,  gravitational  circulation  due  to  freshwater 
discharge,  and  Coriolis  force  due  to  the  earth  rotation.  An  estuary  where  gravita- 
tional circulation  is  dominant  is  called  a salt  wedge  estuary.  But  the  tidal  oscillation 
increases  mixing  between  the  salt  and  freshwater  and  reduces  the  effectiveness  of  the 
gravitational  circulation.  James  River  is  called  a partially  mixed  estuary  because 
both  the  gravitational  circulation  and  the  tidal  mixing  play  an  important  role  in  the 
estuarine  dynamics. 

The  James  River  (Figure  1.1)  originates  in  the  Appalachian  mountains  of  Virginia 
at  the  confluence  of  the  Jackson  and  Cowpasture  Rivers.  The  river  runs  approximately 
530  km  into  Chesapeake  Bay  southeastward.  About  one  third  of  it,  the  160  km  (from 
the  mouth  at  Sewell’s  Point  to  the  fall-line  at  Richmond)  is  a tidal  zone,  and  about 
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one  third  of  this  tidal  zone  is  the  estuary  proper.  The  central  half  of  the  estuarine 
portion  of  the  James  River  contains  some  of  the  most  productive  seed  oyster  beds  in 
the  world.  These  beds  cover  an  area  of  25,000  acres,  of  which  16,150  acres  have  been 
determined  to  be  higher  or  moderately  productive  (Haven  and  Whitcomb,  1983). 

The  average  tidal  range  varies  from  76  cm  at  Sewell’s  Point  to  58  cm  near  the 
Chickahominy  River  mouth  to  98  cm  at  Richmond.  Mean  annual  freshwater  dis- 
charge at  Richmond  is  215  m3/sec  and  in  the  summer  months  flow  is  more  typically 
in  the  range  of  35  to  75  m3/sec.  (Cerco,  1982).  Salt  generally  intrudes  upstream  to 
the  Jamestown  Island,  but  the  intrusion  may  approach  Hopewell  or  be  forced  down- 
stream as  far  as  Newport  News  due  to  low  flow  and  heavy  rainfall.  The  lower  James 
River  near  Hampton  Roads  generally  shows  an  extremely  complex  three-dimensional 
and  time- varying  density  field  due  to  complicated  geometry,  highly  irregular  bottom 
topography,  tidal  flow,  wind  effect,  and  variable  river  discharge. 

As  shown  in  Figure  1.1,  the  James  River  is  channelized  as  far  upstream  as  Rich- 
mond with  a minimum  depth  at  mid-channel  of  approximately  8 m.  Maximum  depth 
in  the  estuary  is  10  - 15  m and  average  depth  is  about  6 m.  Cross-sectional  area 
ranges  from  approximately  27,000  m2  near  the  James  River  Bridge  to  1,000  m2  at 
Richmond.  In  general,  the  basin  is  on  the  order  of  7 km  wide  from  the  mouth  at 
Old  Point  Comfort  to  a constricted  region  at  Mulberry  Point,  just  upstream  of  the 
broader  (10  km  wide)  region  of  Burwell  Bay.  The  remaining  portion  of  the  estuary, 
from  Mulberry  Point  to  Jamestown  Island  is  from  4 to  6 km  wide. 
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figure  1.1:  The  James  River  Estuary  from  Old  Point  Comfort,  at  the  mouth,  to 
Jamestown  Island  (upper  left  corner).  Seed  oyster  beds  are  shown  as  irregularly 
shaped  darkened  areas  (Byrne  et  al.,  1987). 


CHAPTER  2 

GOVERNING  EQUATIONS  AND  SOLUTION  TECHNIQUE 


In  this  chapter,  the  basic  mathematical  model  for  estuarine  circulation  is  de- 
scribed. First  of  all,  assumptions  and  differential  equations  for  estuarine  circula- 
tion models  are  described  in  terms  of  a Cartesian  (a;,  y,  z ) coordinate  system.  The 
equations  are  then  transformed  into  a vertically-stretched  ( x,y,cr ) coordinate  system. 
Following  that,  equations  in  dimensionless  form  are  presented.  Anticipating  model 
applications  to  estuaries  with  complex  geometry,  a method  for  generating  boundary- 
fitted  grid  is  briefly  discussed.  Equations  in  the  boundary-fitted  (£,77,<t)  coordinate 
system  are  then  presented.  Following  a presentation  of  boundary  conditions,  initial 
conditions,  solution  technique,  and  finite-difference  equations  are  presented. 

2.1  Differential  Equations 

The  four  basic  assumptions  used  here  for  deriving  the  governing  equations  are: 
(1)  The  flow  is  incompressible  so  that  the  flow  is  nondivergent.  (2)  The  vertical 
accelerations  are  negligible  compared  to  gravity,  thus,  the  vertical  pressure  distribu- 
tion satisfies  hydrostatic  approximation.  (3)  Density  is  approximated  by  its  mean 
value  except  when  multiplied  by  gravity,  t.e.,  the  Boussinesq  approximation  is  valid. 
(4)  The  concept  of  turbulent  (or  eddy)  viscosity  and  diffusivity  is  introduced.  The 
smallest  scales  of  motion  in  real  fluids,  i.e.,  the  dissipative  scales,  are  always  several 
orders  of  magnitude  larger  than  the  molecular  scales  (Tennekes  and  Lumley,  1972). 
In  contrast  to  their  molecular  parts,  turbulent  viscosity  and  diffusivity  are  not  fluid 
characteristics  but  depend  strongly  on  the  state  of  turbulence.  They  are  properties  of 
the  flow  instead  of  the  fluid  and  hence  are  temporally  and  spatially  dependent.  Us- 
ing a right-handed  Cartesian  coordinate  system  (x,y,z)  with  the  origin  at  the  water 
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surface  and  z measured  upward  as  shown  in  Figure  2.1,  the  equations  may  be  written 
as  follows  (Sheng,  1989) 


du  dv  dw 
dx  ^ dy  ^ dz 


(2.1) 


du  du 2 duv  duw  1 dp  d ( du\ 

dt  dx  + dy  + dz  V p0dx  dx  \ H dx ) 

d (a  du  \ d ( A du' 

+ dy  ( H dy  ) + dz  ( v dz 


(2.2) 


dv  duv  dv2  dvw 
dt  ^ dx  ^ dy  ^ dz 


l dp  d(  dv' 
p0dy  + dx  \ H dx, 
d ( . dv\  d ( . dv\ 

+ a?  r"5i/) + 5 Ir&) 


(2.3) 


dp 

Tz  = ~P9 


dT  duT  dvT_  dwT 
dt  ^ dx  ^ dy  ^ dz 


d ( dT\  d ( dT' 

dx  \ H dx  ) dy  \ H dy , 

df  dT\ 

+ dz  ( Kvdz) 


(2.4) 

(2.5) 


dS  duS  dvS  dwS 
dt  ^ dx  ^ dy  ^ dz 


d / dS\  d / dS' 

dx  y H dx  ) dy  \ H dy , 

d (n  dS\ 

+ dz[Dvdz) 


(2.6) 


p = p(T,S)  (2.7) 

where  (u,v,w)  are  mean  fluid  velocities  in  the  (x,y,z)  directions,  / is  the  Coriolis 
parameter  defined  as  2 fi  sin  (f>  where  fi  is  the  rotational  speed  of  the  earth  and  <f>  is 
the  latitude,  g = 9.81ms-2  is  the  Earth’s  gravitational  acceleration,  p is  density,  p is 
pressure,  T is  temperature,  S is  salinity,  (AH,  KH,  DH)  are  horizontal  turbulent  eddy 
coefficients,  and  (Av,  Kv,  Dv)  are  vertical  turbulent  eddy  coefficients.  Variable  eddy 
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Nominal  Water 
Surface (z=0) 


Figure  2.1:  The  right-handed  Cartesian  coordinate  system. 


Temperature  (°C) 
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Salinity  (°/oo) 


Figure  2.2:  at  density  field  calculated  from  Eckart  formula. 
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coefficients  are  allowed  in  the  present  model.  In  addition  to  the  above  equations,  the 
program  code  contains  an  equation  for  dissolved  species  concentration  similar  to  Eq. 
(2.6). 

Equation  (2.1)  represents  the  continuity  equation  as  conservation  of  water  mass. 
Equations  (2.2)  and  (2.3)  are  conservation  of  momentum  in  the  x and  y directions. 
Eq.  (2.4)  expresses  the  hydrostatic  balance.  Equations  (2.5)  and  (2.6)  represent  a 
conservation  of  energy  and  conservation  of  salt,  respectively.  An  equation  of  state  for 
seawater,  Eq.  (2.7),  relating  salinity  and  temperature  to  density  closes  the  system. 
In  the  present  model,  the  equations  given  by  Eckart  (1958)  are  used 

p = (P  + l)/(a  + 0.698P) 

P = 5890  + 38P  — 0.375T2  -f  3S  (2.8) 

a = 1779.5  + 11.25T  — 0.0745T2  — (3.8  + 0.0171)5' 

where  T is  in  °C,  S is  in  %o,  and  p is  in  g cm-3.  Figure  2.2  shows  crt  density  field 
calcualted  from  Eckart  formula. 

2.2  Governing  Equations  in  Stretched  Coordinates 

To  simulate  circulation  and  transport  in  water  bodies  with  gradual  bathymetric 
variation,  it  is  common  to  use  cr-transformation  such  that  both  the  free  surface  and  the 
bottom  become  the  coordinate  surfaces  with  an  equal  number  of  coordinate  surfaces 
in  between.  This  transformation,  be.,  the  so  called  cr-stretching,  leads  to  a smooth 
representation  of  the  topography  with  the  same  order  of  vertical  resolution  for  the 
shallow  and  deeper  parts  of  the  water  body  (Figure  2.3).  The  cr-stretching  introduces 
extra  terms  to  the  original  Cartesian  equations  of  motion.  However,  most  of  these 
extra  terms  appear  in  the  horizontal  diffusion  terms,  which  are  generally  not  very  im- 
portant compared  to  the  other  terms.  The  cr-grid  model  is  suitable  for  simulating  flow 
and  salinity  transport  in  regions  of  gradual  bathymetric  variation  and  gives  a more 
accurate  estimation  of  bottom  stress  than  the  z-grid  model,  which  resolves  the  depth 


with  “stair-step”  grids.  Recent  studies  by  Sheng  et  al.  (1989c)  and  Haney  (1990) 
found  that,  when  there  are  sufficient  grid  points  across  regions  of  sharp  bathymetric 
change,  the  cr-grid  model  can  yield  good  results.  In  case  of  insufficient  grid  points, 
special  care  must  be  exercised.  Sheng  et  al.  (1989c)  suggested  direct  evaluation  of 
the  horizontal  density  gradient  terms  along  the  constant  2-plane  and  avoiding  the  use 
of  higher-order  advective  schemes  along  the  sharp  bathymetric  variation  to  reduce 
numerical  error.  This  issue  will  be  further  discussed  later. 

The  transformation  from  the  Cartesian  grid  (x,y,z,  t)  to  the  vertically-stretched 
grid  (x7,  y7,  a,  t')  is  defined  as 


/ / ( i\  z C(a'5  2/?  ^)  J j 

x = x,  y = y,  a(x,y,z,t ) = ■■  ■■  ■ ,t  = t 


(2.9) 


h(x,y)  + C{x,y,t)’ 

where  h is  the  water  depth  relative  to  mean  sea  level,  and  ( is  the  free  surface  ele- 
vation relative  to  mean  sea  level,  and  a is  the  transformed  vertical  coordinate  such 
that  a = 0 at  the  free  surface  and  a = 1 at  the  bottom.  Using  the  chain  rule,  the 
following  partial  differential  expression  can  be  derived 


dF__dF_dx^_  d£da  dF  dt' 

dx  dx'  dx  dy'  dx  da  dx  dt'  dx 


(2.10) 


where  F(x , y,  z , t ) = F(x',  t/',  a,  t').  The  first-order  differential  operators  in  the  (x,  y,  z , t ) 
coordinate  system  are  related  to  those  in  the  (x7,  y7,  cr,  t ) coordinate  system  as 


d 

d 

d 

dx 

dx' 

lxda 

d 

d 

d 

dy 

dy' 

lyda 

d 

1 d 

dz 

H da 

d 

d 

1 ,/ 

dt 

dt'  ~ 

(2.11) 


df'ila 


where 
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Figure  2.3:  The  vertically  stretched  grid  (From  Sheng  and  Butler,  1982). 
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The  hydrostatic  equation  in  the  transformed  coordinate  system  is 


(2.12) 


dp  = —pgdz  = —pgHda 
Integrating  from  the  surface  (a  = 0)  to  some  depth  <r, 


(2.13) 


=»>.■-»/  pH  da  (2.14) 

Jo 

where  pa  is  atmospheric  pressure.  Taking  the  ^-derivative  while  applying  the  coordi- 
nate transformation  results  in 


r\££__LfKj.  , oq  fom 

dx  dx'  9 Jo^dx'  H^dx'  + a dx'  do  a + 9p0dx'  ^2'15^ 

where  p0  is  the  mean  density.  After  some  algebraic  manipulation,  Eq.  (2.15)  can  be 
written  as 


dp  l.d(  dH  dp 


dp  _ dpa  rdpH dH  , d( 
'Wd‘T+9d?'7p  + 9d* 


(2.16) 


Now  we  will  drop  all  the  primes^)  for  convenience.  Using  the  flux  conservative 
form  for  the  nonlinear  advective  terms  and  neglecting  the  atmospheric  pressure  gra- 
dients, complete  equations  of  motion  can  be  written  as  follows  (higher-  order-terms 
(H.O.T.)  are  shown  in  Appendix  A) 

Conservation  of  momentum  in  x direction 


dH 

dt 


+ 


u dHu 2 dHuv  dHuu  H r0 

~dr+~^r+~d^  = fvH+!'7jJo 


- »*£+< 


dpH  , dH . 
do  — pa -^— ) 


dx 

1 d 


dx 


(2.17) 


du 


+ Jp  )]  + H.O.T. 


da 
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Conservation  of  momentum  in  y direction 


dHv 

dt 


+ 


dHuv  dHv 2 dHviv  TI  H r dpH 


dx 


+ 


dH , 


= -fuH  + g—(  [ da  - pa^)  (2.18) 
Po  Jo  dy  dy 

TId(  d ( dv\  d ( dv\  1 d ( dv\ 
gHd^  + H^  +dy  \Ah dy)  + ~H*fo  \Av Ik) ' '* 


Conservation  of  water  mass 


(K  dHu  dHv^  dHu 
dt  ^ dx  ^ dy  ^ da 


Conservation  of  heat 


dHT  dHuT  dHvT 

+ —z — + —z — + 


dt 


dx 


dy 


Conservation  of  salt 


dt 


dx 


dy 


where  the  vertical  velocity,  u,  is 


dHuT  d ( dT' 
da  ~H  dx[  <H  dx  , 


TT  d / dT\  1 d / dT' 
+ Hdy\hHdy)  + Hda\Kvda/ 


dHS  dHuS  dHvS 
+ — — + — — + 


dHuS  ud  /_  dS s 
da  dx  1 H dx  , 


„ d / dS\  1 d / dS' 
+ dy{  Hdy)  + Hda{°vda, 


(2.19) 


(2.20) 


(2.21) 


da  w 1 


A 


U=di=H-H{l+a)Tt-^U-1‘V 

which  is  a result  of  the  a coordinate  transformation. 


(2.22) 


2.3  Dimensionless  Equations  in  Stretched  Coordinates 


The  dimensionless  equations  make  it  easy  to  compare  the  relative  importance  of 
various  terms  in  the  equations.  In  a properly  non-dimensionalized  equation  to  be 
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shown  later,  the  part  of  each  term  contained  within  a parenthesis  or  bracket  should 
be  of  order  unity,  while  the  part  of  that  term  containing  the  dimensionless  number(s) 
should  indicate  the  order  of  the  term.  This  should  allow  one  to  compare  the  rela- 
tive importance  of  all  the  terms  and  to  simplify  the  model  equations  under  certain 
circumstances.  The  following  non-dimensionalization  procedure  follows  that  used  by 
Sheng  (1986b)  and  Sheng  et  al.  (1991).  The  added  insight  justifies  the  slightly  more 
complicated  form  of  the  dimensionless  equations. 

Dimensionless  Variables 


(u*,v*,w*) 

= 

(U,V,WXTIZT)  /UT 

(x*,y*,z*) 

= 

(x,y,zXr/Zr)  /Xr 

(TriTy) 

= 

(Tx  5 Ty)/ {PofZrUT)  = (r*,r; 

t* 

= 

tf 

Q*s 

— 

T0/(Tr  - T0)  * q/(p0cpfZrT0 ) 

c 

— 

gC/(fUrXr)  = C/ Sr 

p* 

= 

(P~  Po)/(pr  ~ Po) 

T* 

= 

(T  — T0)/(Tr  — T0) 

A*h 

= 

a: 

= 

Ay  / Ayr 

K'h 

= 

KH/KHr 

K 

= 

Kv/Kyr 

D'„ 

=7 

DhI  DHt 

d: 

= 

Dy/Dyr 

U)* 

= 

uXr/Ur 

(2.23) 
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Dimensionless  Parameters 


Vertical  Ekman  Number: 
Lateral  Ekman  Number: 
Vertical  Prandtl  Number: 
Lateral  Prandtl  Number: 
Vertical  Schmidt  Number: 
Lateral  Schmidt  Number: 
Froude  Number: 
Rossby  Number: 


Ey  = Avr/(fZ? ) = ti/tvdm 


Eh  = Aht/UX?)  = ti/tidm 

Pf  v — Avr  / Kvr  — tvdh  / tydm 
Prn  — A Hr  I K Hr  = tedh/tedm 


Scv  — Avr/Dvr  — tvds/lvdm 


ScH  = AHt/ DHr  = tldjtldm  (2.24) 


Fr  = UT/(gZTy /2  = tge/tc 
Ro  = Urf{fXr)  = ti/tc 


P = gZr/(pX r2)  = ( Ro/Frf  = (U/tgey 
Sz  = Sr/Zr 

£ = (Pt  — Po)/po 

Densimetric  Froude  Number:  Fro  = Fr/y/e  = tgi/tc 


where  the  various  V s appearing  in  the  above  expressions  represent  the  characteristic 
time  scales  associated  with  various  physical  processes  in  lakes  and  estuaries  (Ta- 
ble 2.1). 

Utilizing  the  dimensionless  variables  and  parameters  defined  above  and  suppress- 
ing the  asterisks  for  clarity,  one  can  derive  the  following  dimensionless  equations 


dC,  adHu  adHv  nrrdu 
4 + P-fi-  + P-7T-  + PH—  = 0 
dt  dx  ay  ocr 


(2.25) 


1 dHu 
H dt 


d(  Kd_  / dvP 
dx  + H2da  \ vda/ 
Ro  ( dHuu  dHuv 
~H  { dx  + dy 


| + v 
dHuu 

+ 


(2.26) 


da 
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Table  2.1:  Characteristic  Time  Scales  of  Various  Physical 
Processes  in  Estuaries  and  Lakes. 


Physical 

Process 

Time 

Scale 

Order  of 
Magnitude 

Periodic  Forcing 

lf 

1 /w 

Convection 

to 

Xr/U, 

Inertia  Oscillation 

u 

1// 

Vertical  Turbulent 
Diffusion 

t vdm ) t vdh  1 1 vds 

Z?/A„,Z2JK„,Z2r/D„ 

Lateral  Turbulent 
Diffusion 

ttdm  i t^dhi  lids 

X?/AHr,X2JKHr,X2rIDHr 

Gravity  Wave 

tge 

Xr/Vez: 

Internal 
Gravity  Wave 

tgi 

Xr  / \J q Zr^p  I p0 

28 


1 dHv 
H~df 


1 dHT 

77  dt 


1 dHS 
H dt 


+ Eh 

Ro 


d (a  9u 
dx  \ H dx 


9 ( A 9u ' 

+ dy  v 


+ //.o.r. 


Erl 


TJ  [°  dp  . dH  ( r°  \ 

hL  rxda + ^ [L pda + ap) 


_ 9C_.E}L9_(a  du\,R 

dx  + H2  da  [Av da  ) + x 


d(  ^Ev  d ( dv\ 

dy  + H2da\Avda)  U 

Ro  ( dHuv  dHvv  dHvu\ 

H ^ dx  dy  da  J 

+ E«\jr  (A«7r)  + T + H 0 T- 

I dx  \ dx  I dy  V dy  J 


Ro 

M>  L 


TT  dp  . dH  ( r°  \ 

H l (L  pd‘T  + V 

d»\By 


_ * . Ey  d 

dy  H2  da  v \da 


(2.27) 


Ev  9 K ?T\ 

PrvH2da\  v da  ) 

Ro  (dHuT  dHvT  dHuT\ 
H \ dx  dy  da  J 

Eh 


Pr 


H 


Ukh^ 

dx  \ dx 


+Uk4)+h-°-t- 


(2.28) 


(D  dA) 

ScvH2da\  v da  ) 

Ro  SdHuS  dHvS  dHuS\ 
H \ dx  dy  A da  ) 

Eh 


Sch 


l < 


+ £ c»f  + HOT- 


(2.29) 


P = p(T,  S ) 


(2.30) 


where  H = h(x,y)  + ((x,y,t)  is  the  total  depth.  The  H.O.T.s  are  omitted  here  for 
simplicity. 
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2.4  Dimensionless  Vertically  Integrated  Equations 

The  external  mode  consists  of  the  surface  elevation  and  the  vertically-integrated 
velocities.  The  vertically  integrated  equations  represent  the  “external  mode”  of  the 
3-D  model.  The  dimensionless  vertically  integrated  equations  are 


dC,  a(dU  dV\  n 


(2.31) 


dU 

dt 


HdC  _ T/ 

" n T T sx  Tbx  T V 

OX 


Ro 


d 

(UU\ 

, d 

(UV\  1 

dx 

UJ 

+ d~y 

H )] 

+ Eh 


d 


dU\  d 


rM»«r  + 


dx 


dx  j dy 


' dU' 


Ro  H2  dp 

FrD2^Tdi 

H^  + Dx 

ox 


(2.32) 


dV 

dt 


d( 

T Tsy  — Tky  — U 


Ro 


d (UV\  d (W\ 

~dx  lirj  + dy  [if) 


d_ 

dy 


'VV' 


+ E1 


H 


d 


dv 


d 


o-  ^ [A 


dV ' 


dx 


dx 


dy  \ H dy 


Ro  H2  dp 
FrD2  2 


(2.33) 


where  ({/,  V)  = /° x H(u,v)dcr  are  the  vertically  integrated  velocities.  The  nonlinear 
inertia,  lateral  diffusion,  and  baroclinic  pressure  gradient  terms  in  Equations  (2.32) 
and  (2.33)  are  obtained  by  vertically  integrating  the  corresponding  terms  in  Equations 
(2.26)  and  (2.27),  respectively.  However,  these  terms  contain  idealized  assumptions 
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about  the  vertical  profiles  of  horizontal  velocities  and  density,  i.e.,  the  horizontal 
velocities  and  density  are  assumed  to  be  uniform  in  the  vertical  direction,  and  hence 
are  not  exact. 

2.5  Boundary-Fitted  Coordinate  Generation 

The  basic  idea  of  the  boundary-fitted  coordinate  system  approach  is  to  make 
computations  on  curvilinear  grids  so  that  one  of  the  curvilinear  coordinates  always 
follows  a boundary.  According  to  the  work  of  Thompson  et  al.  (1974),  the  natural 
coordinates  £ and  rj  are  taken  as  solutions  of  an  elliptic  boundary  value  problem  with 
one  of  the  coordinates  being  constant  on  the  boundaries.  The  same  methodology  was 
used  to  generate  the  curvilinear  grids  for  the  present  study.  Detail  discussions  are 
given  Thompson  et  al.  (1974,  1976,  1977a,  b),  Thames  (1975),  and  Thames  et  al. 
(1977). 

The  curvilinear  coordinates  are  determined  by  solving  an  elliptic  system 

(xx  + Zyy  = Pti,r))  (2.34) 

IJxx  T Vyy  — Q(£iV)  (2.35) 

with  Dirichlet  boundary  conditions, i.e.,  the  values  of  £ and  r]  are  specified  as  constants 
on  the  boundaries.  ( x,y ) is  coordinate  system  in  the  physical  domain,  and  (£,t?)  is 
coordinate  system  in  the  transformed  domain.  Subscripts  indicate  differentiation.  To 
obtain  the  grid  locations  in  the  (x,y)  plane  corresponding  to  known  uniform  rectan- 
gular grid  in  the  (£,  rj)  plane,  the  dependent  and  independent  variables  in  Eq.  (2.34) 
and  Eq.  (2.35)  are  interchanged.  This  results  in  a coupled  system  of  quasi-linear 
elliptic  equations  for  determining  a:(£,7/)  and  y(£,rj)  in  the  transformed  plane.  Fig- 
ure 2.4  shows  the  relation  between  physical  domain  and  transformed  domain. 


axU  - Wxlr,  + lxr1ri  + aX^P  + ^XVQ  = 0 


(2.36) 
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Physical  Plane 


Transformed  Plane 


Figure  2.4:  Field  transfomation  (From  Thompson  et  al. , 1977). 
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(2.37) 


ay«  - 2 pytr,  + IVriv  + ay^P  + iynQ  = 0 
where  the  coefficients  are 

2 I 2 

« = *;  + y„ 

f3  = X^XTJ  + yiVr) 

7 = ®|  + y? 

Q((,v)  = j2iQ 

J = Xtfr,  - x„y^ 


(2.38) 


Here  P and  Q are  the  control  functions  that  serve  to  concentrate  the  grid  lines. 
With  a weight  function  w,  P = w^/w,  Q = wv/w. 

2.6  Transformation  of  Governing  Equations  Into  Boundary- Fitted  Grids 

2.6.1  Transformation  Rules 

Generation  of  a boundary-fitted  grid  is  an  essential  step  in  the  development  of 
a boundary-fitted  hydrodynamic  model.  It  is,  however,  only  the  first  step.  A more 
important  step  is  the  transformation  of  governing  equations  into  the  boundary-fitted 
coordinates.  A straightforward  method  is  to  only  transform  the  dependent  variables, 
i.e.,  the  coordinates,  while  retaining  the  Cartesian  components  of  velocities.  John- 
son (1982)  developed  such  a two-dimensional  vertically-integrated  model  of  estuarine 
hydrodynamics.  The  advantage  of  such  method  is  its  simplicity  in  generating  the 
transformed  equations  via  the  chain  rule.  The  dimensional  forms  of  the  continuity 
equation  and  the  vertically-  integrated  momentum  equations  are  shown  by  Equations 
(2.17),  (2.18)  and  (2.19).  The  resulting  equations,  however,  are  rather  complex.  Even 
when  an  orthogonal  or  a conformal  grid  is  used,  the  equations  do  not  become  any 
simpler.  Additional  disadvantages  are:  (1)  the  boundary  conditions  are  quite  compli- 
cated because  the  Cartesian  velocity  components  are  generally  not  aligned  with  the 
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grid  lines,  (2)  the  staggered  grid  cannot  be  readily  used,  and  (3)  numerical  instability 
may  develop  unless  additional  variables  (e.g.,  surface  elevation  or  pressure)  are  solved 
at  additional  grid  points  (Sheng,  1989). 

To  alleviate  the  problems  mentioned  in  the  previous  paragraph,  we  have  chosen 
to  transform  the  dependent  variables  as  well  as  the  independent  variables.  Equations 
in  the  transformed  coordinates  (£,r))  can  be  obtained  in  terms  of  the  contravariant, 
or  covariant,  or  physical  velocity  components  via  tensor  transformation  (Sokolnikoff, 
1960).  Rules  of  transformation  and  dimensionless  tensor  invariant  governing  equations 
are  described  in  Appendix  B. 

2.6.2  Dimensionless  Equations  in  Boundary-Fitted  Grids 


It  is  clear  that  if  a proper  transformation  is  given  by  Equations  (2.36),  (2.37), 
(2.38)  and  (2.39),  one  can  compute  the  metric  tensor  and  Christolfel  symbols  accord- 
ingly, and  then  expand  Equations  (2.48)  and  (2.49)  in  terms  of  the  two  contravariant 
velocity  components  and  sum  up  all  the  terms  over  the  range  of  indices.  This  process 
is  rather  tedious  and  highly  prone  to  human  error.  To  alleviate  these  problems,  we 
used  a symbolic  manipulator  to  expand  the  tensor-invariant  equations.  After  rear- 
ranging the  resulting  terms  of  expansion,  we  obtain 
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+ £//A//(Horizontal  Diffusion) 
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where  the  horizontal  diffusion  terms  are  listed  in  Appendix  C.  The  temperature  and 
salinity  equations  can  be  obtained  according  to  the  same  procedure  as 
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Recently,  the  nonlinear  terms  in  the  momentum  equations  have  been  replaced  by  a 
more  compact  formulation  (Sheng,  1989).  The  u-equation  and  u-equation  are  shown 
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in  Appendix  D.  The  more  compact  equations  do  yield  better  numerical  stability  than 


the  earlier  equations  shown  in  Equations  (2.40)  and  (2.41). 

2.7  Boundary  Conditions  and  Initial  Conditions 


Boundary  conditions  are  required  for  computing  the  six  dependent  variables  on 


or  near  the  boundary.  If  the  closed  boundary  is  fixed  in  space,  the  component  of 


to  the  boundary).  If  the  fluid  is  viscous,  an  additional  condition,  i.e.,  the  no-slip 


the  closed  boundary,  no  normal  flux  of  momentum,  salt,  or  energy  is  allowed  since 
the  velocity  is  zero.  All  the  variables  shown  below  are  dimensionless  based  on  the 
nondimensionalization  in  section  2.3  except  for  S. 

2.7.1  Vertically  Stretched  System 

The  boundary  conditions  at  the  free  surface  (a  = 0)  are 


fluid  velocity  normal  to  it  must  be  zero  (V  ■ n = 0,  where  n is  the  unit  normal  vector 


condition  (V- s = 0,  where  sis  the  unit  tangential  vector  to  the  boundary)  is  applied. 
Application  of  the  above  closed  and  no-slip  conditions  is  made  here  which  result  in 
the  fluid  at  the  boundary  having  zero  velocity  with  respect  to  the  boundary  itself.  At 
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The  boundary  conditions  at  the  bottom  (cr  = —1)  are 


(2.45) 


(2.46) 
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i-°  <2-47> 

where  u\  and  Vi  represent  the  velocity  components  at  the  first  grid  point  above  the 
bottom. 

Along  the  shoreline  where  river  inflow  or  outflow  may  occur,  the  conditions  are 
generally 


u = u(x,y,cr,t) 
v = v(x,y,a,t ) 

u>  = 0 (2.48) 

T = T(x,y,<r,t) 

S = S(x,y,a,t) 


The  kinematic  boundary  condition  is  also  imposed  at  the  surface. 


d(  d(  d( 
u,=  m+uai  + v9i 


(2.49) 


The  boundary  condition  required  for  the  pressure  is  P = 0 at  <r  = 0 which 
expresses  the  assumption  that  either  the  contribution  of  the  atmospheric  pressure 
to  the  pressure  of  the  fluid  column  is  negligible,  or  that  the  atmospheric  pressure 
gradients  are  unimportant. 

The  open  boundary  condition  is  much  more  complicated,  since  momentum  and 
mass  transport  occur  through  the  boundary.  Salinity  and  temperature  need  to  be 
specified  with  £ or  velocity  at  the  open  boundary.  For  the  elevation  £,  there  are 
currently  several  options 
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where  An,  Tn  and  <j)n  are  the  amplitude,  period  and  phase  angle  of  the  tabular  tidal 
constituents.  When  open  boundary  conditions  are  given  in  terms  of  (,  the  normal 
velocity  component  is  assumed  to  be  of  zero  slope  while  the  tangential  velocity  com- 
ponent may  be  either  (1)  zero,  or  (2)  zero  slope,  or  (3)  computed  from  the  momentum 
equations.  But  in  reality  the  velocity  gradients  are  impossible  to  obtain  from  field 
observation.  In  the  present  study,  it  is  assumed  that  the  velocity  gradients  are  zero 
at  the  open  boundary. 

The  salinity  along  an  open  boundary  or  river  entrance  is  generally  prescribed 
throughout  the  tidal  cycle.  The  present  model  allows  the  salinity  at  the  open  bound- 
ary to  be  computed  from  a 1-D  advection  equation  during  the  outflow.  For  example, 
along  an  open  boundary  on  the  left  or  right  (Sheng,  1986) 

dHS  „ dHuS 

-dT  + R°-ar 

where  the  spatial  salinity  flux  is  evaluated  from  the  salinity  values  at  the  boundary 
and  the  interior  grid  point  via  a one-sided  differencing  scheme.  During  the  inflow, 
however,  the  salinity  value  at  the  open  boundary  can  either  take  on  a prescribed  value 
or  be  determined  from  the  1-D  advection  equation  while  using  the  boundary  salinity 
value  and  the  prescribed  salinity  value  to  evaluate  the  spatial  flux  term. 

To  initiate  a run,  the  initial  spatial  distributions  of  £,  u,  v , oj,  T and  S need  to 
be  specified.  When  initial  data  are  unknown,  one  could  start  with  zero  initial  fields. 
When  initial  data  are  known  at  a limited  number  of  locations,  an  initial  field  can 
be  generated  by  some  interpolation  scheme.  It  is  desirable  if  the  interpolated  field 
satisfies  the  conservation  law  governing  that  field  variable. 

For  practical  simulations,  the  present  model  usually  assumes  zero  initial  flow  if 
little  initial  data  are  known.  In  case  of  salinity  simulation,  the  spin-up  time  is  longer 
and  an  interpolation  routine  is  provided  to  produce  a reasonable  initial  field  from 
limited  data  points. 
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2.7.2  Vertically  Stretched  and  Laterally  Boundary-Fitted  System 
The  boundary  conditions  at  the  free  surface  (a  = 0)  are 
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The  boundary  conditions  at  the  bottom  (a  = — 1)  are 
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where  Uj  and  iq  are  the  contravariant  velocity  components  at  the  first  grid  point 
above  the  bottom. 

Due  to  the  use  of  contravariant  velocity  components,  the  lateral  boundary  condi- 
tions in  the  (£,rj,a)  grid  are  similar  to  those  in  the  (x,y,a)  system.  Along  the  solid 
boundary,  the  no-slip  condition  dictates  that  the  tangential  velocity  is  zero  while  the 
slip  condition  requires  that  the  normal  velocity  is  zero.  When  flow  is  specified  at  the 
open  boundary  or  river  boundary,  the  normal  velocity  component  is  prescribed. 

Initial  conditions  on  vectors,  if  given  in  the  Cartesian  or  prototype  system,  such 
as  the  velocity  and  the  surface  stress,  must  be  first  transformed  before  being  used 
in  the  transformed  equations.  Thus,  the  surface  stress  in  the  transformed  coordinate 
system  is  given  by 
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where  r1,  r2  are  the  contravariant  components  of  the  stress  in  the  transformed  system 
and  r1,  r2  are  the  contravariant  components  in  the  Cartesian  system.  Note  that 
in  the  Cartesian  system,  the  contravariant,  covariant  and  physical  components  of  a 
vector  are  identical.  The  contravariant  components  of  the  initial  velocity  vectors  can 
be  transformed  in  the  same  manner  to  obtain  the  proper  initial  conditions  for  the 
transformed  momentum  equations. 

2.8  Solution  Technique 

The  solution  of  the  difference  equations  may  differ  from  the  solution  of  the  differ- 
ential equations  because  of  a variety  of  errors  (truncation  error,  dispersion  error,  and 
round-off  error)  even  if  the  scheme  is  stable.  For  explicit  methods,  the  stability  con- 
dition is  that  the  CFL  (Courant-Friedrichs-Lewy)  criteria  based  on  the  shallow  water 
wave  speed  (Roach,  1976)  should  be  satisfied.  The  basic  numerical  technique  to  en- 
hance the  model  efficiency  is  the  mode-splitting  technique  (Sheng  et  al.,  1978;  Sheng, 
1985;  Madala  and  Piascek,  1977),  which  basically  solves  the  equations  of  motion  by 
splitting  them  into  an  “external”  or  vertically-averaged  mode  and  an  “internal”  or 
vertical  structure  mode.  The  external  mode  is  solved  implicitly,  and  the  internal 
mode  is  solved  explicitly  except  for  the  vertical  diffusion  terms.  Semi-implicit  meth- 
ods (Sheng,  1983;  Sheng  and  Choi,  1989)  can  be  introduced  to  take  advantage  of  the 
merits  of  both  the  explicit  method  and  the  implicit  method.  This  approach  solves 
implicitly  those  terms  in  the  governing  equations  which  possess  rather  stringent  sta- 
bility criteria  if  explicit  schemes  were  used.  These  terms  include  those  that  govern  the 
propagation  of  surface  gravity  waves:  local  acceleration,  and  surface  slope  terms  in 
the  momentum  equations  and  all  terms  in  the  continuity  equation.  Furthermore,  ver- 
tical diffusion  terms  in  the  momentum  and  transport  equations  are  solved  implicitly. 
All  other  terms  in  the  equations  are  solved  explicitly. 

The  governing  equations  are  solved  by  a finite  difference  scheme  along  with  proper 
boundary  conditions.  A space-staggered  grid  system  shown  in  Figure  2.5,  is  used  to 
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discretize  the  differential  equation  (Sheng,  1983).  In  the  staggered  grid,  elevation, 
vertical  velocity,  temperature,  salinity,  and  density  are  defined  at  the  center  of  a cell, 
and  it- velocity  at  the  left  and  right  sides  of  a cell,  and  u-velocity  at  the  top  and  bottom 
sides  of  a cell.  A staggered  grid  system  has  an  accuracy  of  0( Ax2)  with  the  same 
mesh  system  compared  with  an  accuracy  of  0(A;r)  in  a non-staggered  grid  system. 
For  the  time  marching,  forward  time  difference  (two-time-level)  or  leapfrog  technique 
(three-time-  level)  can  be  used.  The  present  study  uses  a two-time-level  scheme. 
Both  the  leapfrog  technique  and  spatial  difference  are  centered  differences,  and  the 
accuracy  in  space  and  time  is  second-order.  The  separation  of  solutions  caused  by 
the  leapfrog  method  is  corrected  by  using  a periodic  time  smoothing.  The  diffusion 
terms  are  evaluated  at  the  previous  time  level  to  avoid  numerical  instabilities. 

In  the  computational  procedure,  first  the  spatial  derivatives  x^^xv,y^yv,  and  Ja- 
cobian J at  each  grid  point  are  computed  using  the  transformation  equations  (Eq. 
(2.34)  - Eq.  (2.38)).  These  derivatives  are  used  as  input  to  the  circulation  and  trans- 
port model.  The  “external  mode”  equations  are  solved  using  an  alternating-direction 
implicit  (ADI)  method.  This  solution  technique  consists  of  a two-step  calculation 
which  advances  the  solution  first  from  time  level  n to  * (an  intermediate  step)  and 
then  from  * to  n + 1.  In  the  first  step,  £ and  U are  advanced  from  n to  * by  an  implicit 
calculation  along  each  row  in  the  £ direction,  while  treating  V and  /^-derivative  terms 
explicitly.  In  the  second  step,  ( and  V are  advanced  from  * to  n + 1 by  an  implicit 
calculation  along  every  column  in  the  rj  direction,  while  treating  U and  ^-derivatives 
explicitly.  Then  deviations  of  the  velocities  from  the  vertically  averaged  velocities  are 
calculated,  while  treating  implicitly  the  vertical  diffusion  terms  in  the  equations  for 
the  deficit  velocities.  The  deficit  velocities  are  then  added  to  the  vertically- averaged 
velocities  to  produce  the  three-dimensional  velocity  field.  Conservation  equations  for 
salinity  and  heat  are  then  solved  explicitly,  except  for  the  vertical  diffusion  term. 
The  implicit  schemes  used  in  the  external  mode  and  the  internal  mode  require  the 
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Figure  2.5:  Staggered  numerical  grids  in  lateral  and  vertical  directions  (From  Sheng, 
1983). 
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inversion  of  only  tridiagonal  matrix  systems. 

2.9  Finite  Difference  Equation 

In  the  following,  a set  of  two-time-level  difference  equations  in  (x,y,a)  system 
are  presented.  A set  of  three- time- level  difference  equations  in  (£,7,  <r)  system  are 
presented  in  Appendix  E.  The  results  presented  in  this  study  are  obtained  with  the 
two- time- level  equations.  The  finite  difference  representation  of  a dependent  variable 
can  be  written 

= <Kx,y,<r,t)  = <f>(iAx,jAy,kAa,nAt) 

6X4>  = (<^+1/2  - ^.-1/2)/ Ax 

Finite-difference  equations  for  the  external  mode  and  the  internal  mode  are  listed 
in  the  following.  The  basic  two-time-level  schemes  are  similar  to  those  presented  in 
Sheng  (1983)  and  Sheng  (1987).  For  three-time-level  schemes,  refer  to  Appendix  E. 
2.9.1  External  Mode  in  (x,  t/,cr)  (Two- Time- Level  Scheme) 

The  external  mode  consists  of  the  surface  displacement  £ and  the  vertically- 
integrated  velocities  U and  V and  is  described  by  Equations  (2.31)  through  (2.33). 
The  finite-difference  schemes  treat  all  of  the  terms  in  the  continuity  equation  implic- 
itly. The  time  derivatives  and  surface  slopes  in  the  momentum  equations  are  generally 
treated  implicitly,  while  the  bottom  stresses  are  computed  explicitly  from  the  latest 
vertical  profiles  of  horizontal  velocities  assuming  the  presence  of  a turbulent  bottom 
boundary  layer.  In  that  case,  one  can  either  specify  a drag  coefficient  or  a rough- 
ness height,  instead  of  the  conventional  Chezy-Manning  coefficients  which  are  used 
in  many  two-  dimensional  models.  When  used  in  two-dimensional  mode,  however, 
the  model  computes  implicitly  the  bottom  stresses  in  the  u and  v equations  with 
with  Chezy-Manning  formulation.  The  dimensionless  finite-difference  equations  in 
the  following  are 


AFn+1  = Fn  + AtDn 


(2.57) 
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where  superscripts  n + 1 and  n denote  the  present  and  previous  time-step  of  integra- 


tion, At  is  the  time-step,  and 
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with  Dx  and  Dy  as  defined  by  Equations  (2.32)  and  (2.33),  and 
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where  Ax  and  Ay  are  grid  spacings  in 

the  x and  y directions,  respectively. 

For  increased  efficiency  of  computation  as  previously  discussed,  the  ADI  method 
factors  the  matrix  A into  the  sum  of  three  matrices:  the  identity  matrix  I,  a matrix 
Ax , and  a matrix  Ay: 
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(2.62) 
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Ay  = 


0 0 Ctyfiy 

0 0 0 

. 7 A 0 0 


(2.63) 


To  solve  the  factored  version  of  Equation  (2.57),  the  ADI  method  proceeds  by  first 
computing  in  the  x-directions  (the  “x-sweep”),  and  then  computing  in  the  ^-direction 
(the  “y-sweep”).  Carrying  out  the  factorization  of  Equation  (2.57),  and  neglecting 
terms  of  order  At2  yields 
x- sweep: 


[/  + AX\F*  = [I  - Ay]Fn  + A tDn 


y- sweep: 


(2.64) 


[/  + Ay]Fn+1  = F*  + AyFn 
where  F*  represents  the  intermediate  solution 

( ^ \ 

F*  = U * 

\V*  ) 

The  more  detailed  equations  for  the  x-sweep  and  y-sweep  are 
x-sweep: 


(2.65) 


(2.66) 


c*  + ax8xU*  = C-  <*y8yVn 
u*  + 7 *8XC  = un  + AtDnx 
V*  = Vn  - 7 y8yC  + A tDl 


(2.67) 

(2.68) 

(2.69) 


The  C*  and  U*  are  two  unknowns  in  the  x-direction.  The  above  equation  can  be 
readily  solved  by  simple  inversion  of  tridiagonal  matrices  along  all  the  grid  rows. 


y-sweep: 
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(2.70) 

(2.71) 


<n+1  + <*y6yVn+i  =C  + «y8yVn 

Vn+l  + 7sACn+1  = ^ + tAC” 

Un+1  = U*  (2.72) 


The  £n+1  and  yn+1  are  two  unknowns  in  the  y-  direction.  Again,  inversion  of 
tridiagonal  matrices  along  all  the  grid  columns  are  performed  to  solve  the  above 
equation  by  substituting  the  equation  for  V*  from  the  x-sweep  into  the  equation  in 
the  y-sweep. 

2.9.2  Internal  Mode  in  (x,  t/,  a ) (Two- Time- Level  Scheme) 


Defining  deficit  velocities  as  u — u — u and  v = v — v,  we  obtain  the  differential 


equations  for  the  internal  mode  by  subtracting  the  vertically-averaged  momentum 


equations  from  the  three-dimensional  momentum  equations  multiplied  by  H (Sheng, 
1986) 


1 dHu 
H ~dt~ 


= Bx 


Dx  Ev  d 


(2.73) 


1 dHv 
H~dT 


_D1,Kd_ 

y H + H2  da 


Avi{v+i\ 


(2.74) 


The  finite  difference  equations  are 


( Hu)n+1 


{Hd)n+l 


= (Hu)n  + At(HBx  - Dx)n 


= (Hv)n  + At(HBy-Dy)n 
'da  v da 


+ §a4^4(^+5)"+i) 


(2.75) 


(2.76) 


The  vertical  diffusion  terms  in  the  momentum  equations  are  treated  implicitly  to 
ensure  numerical  stability  in  shallow  waters.  It  is  also  important  to  ensure  that  the 
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vertically-integrated  perturbation  velocities  are  always  equal  to  zero 


kmax 


X = 0 

k=l 

(2.77) 
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(2.78) 

In  order  to  ensure  that  Equations  (2.77)  and  (2.78)  are  satisfied,  the  nonlinear 
inertia  and  turbulent  diffusion  terms  in  the  vertically-integrated  Equations  (2.32)  and 
(33)  must  be  evaluated  by  summing  the  corresponding  terms  in  the  three-dimensional 
equations  at  all  vertical  levels. 

Once  un+1  and  un+1  are  obtained,  u and  v can  be  obtained  from 


uB+1  = un+l  +Un+1/Hn+1 
vn+1  = vn+1  + Vn+1/Hn+1 


(2.79) 

(2.80) 


Following  these,  the  vertical  velocities  u>n+1  and  wn+1  can  be  computed  from 
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(2.82) 


The  vertical  velocity  u>  should  be  nearly  zero  at  the  free  surface. 
2.9.3  External  Mode  in  (£,  7 /,cr)  (Two-Time-Level  Scheme) 


The  present  study  generally  uses  the  following  two- time- level  scheme  instead  of  the 
two- time- level  scheme  described  earlier.  Vertically-integrated  differential  equations 
are  written  as  follows  (Sheng,  1986) 


Ut  + Hgu(t  + Hgl\v  = ^=JU  - CTiU  + D'x 

y/9o 
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where  CT^  is  drag  coefficient  in  ^-direction,  and  Dx  and  Dy  are  all  the  terms  in  the 
^-direction  and  y-direction  which  are  not  shown  in  the  momentum  equations.  The 
above  equations  are  written  in  finite  difference  forms  as  follows 


£n+x  _ £n  _j_  >crm«+i 


+ 


A i-y/g~0 
(3AT 
Ay  • sfgl 


[TMVfoU)n+l  + (1  - TtMy/gJT)]  (2.86) 

[TMV9~oVn+x)  + (1  - Ti)Sv(y/g^Vn)\  = 0 


U "+1  - un  + 


+ 


Ar{^  [r, <£«"+')  + (1  - TiMC)]  (2.87) 


Hg 


12 


Ay 
- - AT 


TrSr'iC*1)  + (l  - T^iOl 


T2CT(Un+x  + (1  - T2)CTiUn  - T3-^Un+1 

\j9o 


- (1  -T3)^Un 

y/fjo 


+ Dxn 


vn+1-vn  + at|^ 


12 


W£(C"+1)  + (i-r.)«f(C”)] 


Hg 


+ 


= - A T 


22  r 

- [r,Mcn+1)  + (i  - TiK(C) 


r2cr,i/n+1  + (i  - t2)ct„v"  + n-^v+1 

yj9o 


+ (i-r3)-^Kn 

y/do 


+ d; 


(2.88) 


Instead  of  solving  the  above  equations  directly,  however,  they  can  be  written  in 
ADI  scheme  as  follows 
£-  sweep: 


C + T,a(S((^-JJ-)  = C-Ml-TiMy/TM")  (2.89) 


“A(n/9^”) 
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rj-  sweep: 


TiVs~olu6t(C)  + 


1 -f  A t I T2Ch  — T3 


g 12 

y/fo). 


{y/giu*) 


Txy/gtfW)  + (VfoUn) 

(1  - T^J'S^C)  ~ (1  ~ T1)x/£7%(Cn) 


- (VfoUn) 


AT(1  - T2)CTi  - A<(1  - T3) 


g 12 

s/go. 


+ d: 


(2.90) 


C“+1  + T^S^V"*')  = C + TtaMs/glV") 


(2.91) 


r.v'M!!«Cn+1)  + 


1 + A<  ( T^CVn  + ^3 
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s/gl) 


ls/g~oVn+1 ] (2.92) 


ri^712^(0  - (1  - - (1  - TOv^^O 


+ (v^n-(v^n) 


At(l-T2)Crv  + At(l-T3) 


g 12 


V^. 


+AT 


Tiy/g-ol12Sr,(Cn+1)  + (y/^oUn+1)  = (JgJT)  + T,  v^712<UCn_1 ) (2-93) 

where  7\,  T2  and  T3  are  all  constants  between  0 and  1.  For  example,  the  wave 
propagation  terms  are  treated  explicitly  if  T\  = 0,  but  implicitly  if  T\  = 1.  m can  be 
n — 1 or  n.  Additional  parameters  appearing  in  the  equations  are 


pAt 

ai 

— 

goA( 

PA  t 

av 

— 

5oA7 

11 

HAtg11 

7 

— 

A£ 

12 

HAtg 12 
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— 

A£ 

22 

HAtg 22 

7 

Ag 

- 21 

HAtgu 

7 

= 

Ag 

(2.94) 
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The  system  of  equations  can  be  rewritten  in  the  following  matrix  form 


[A}Fn+1  = AFn~x  + D 

(2.95) 

[I  + At\F*  = [B-Av]Fn-1  + D 

(2.96) 

[I  + Av]Fn+1  = F*  + AvFn~1 

(2.97) 

where 


1 <xxT\Sz  atyTrf,, 

[A]  = f 1 + AT  (Ct(T2  - ^r3)  0 

vW2^,  0 1 + AT  (c„t2  + ^r3) 


(2.98) 


W = 


0 ctxT\8^  0 

v^57 "TiS(  AT  (Ct(T2  - *±Ts)  0 


0 


(2.99) 


[A,]  = 


0 

0 


0 

0 


ayTi8v 

0 


v5o722r,«„  o at(ct,t2  + *lt3) 


(2.100) 


[B}  = 


1 0 0 
0 1 0 
0 0 1 


(2.101) 


c 

F=  y/%U 
. V9~oV 


(2.102) 
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D = 


-0,(1  ~ TMy/gSU*))  - ay(  1 - TMJ&V*) 

y/%UnAT  (-(1  - T2)CTi  + (1  - T3)^=)  - v^7n(l  - T^C 
-y/%i216vC  + A TD'xy/g~0 


-V^7%Cn  + A TD'xy/g~Q 

y/%VnAT  (-(1  - T2)Ctv  - (1  - T3)^=)  - v^712^Cn 
L -x/^722(1  - WC"  + A TD'y^ 

2.9.4  Internal  Mode  in  (£,//,  (7)  (Two-Time- Level  Scheme) 


(2.103) 


JU  a = 

ot  y/gl 


12  22 


Ev_d_  / ^ > 

+ H2  5(7  V ” 5(7  ) 


\[9o 

dHu\ 


(M  ~ Tbi)  ~ F2 


ltm  = 


9n  «21 

:i/u-^HG3 


y/9o 


y/9o 

Ev  5 / 5//u\  . „ 

T rro  o I -'d«  o I (Tmj  ^(ir;)  @2 


(2.104) 


(2.105) 


H2  5(7  V * 5(7  # 

and  F2  and  G2  indicate  all  the  explicit  terms  in  U and  V equations,  respectively,  and 
F3  and  G3  indicate  all  the  explicit  terms  in  the  u and  v equations,  respectively. 


,12 


(Hu)n+1  = (Hu)n  + 2At—Hn+1un+1 

y/9o 

q22 

+ 2At- — Hnun  + 2At(F3  - F2)n 

yJ9o 


+ 2At{Hn)28a 

~ (T*<  ~ Tb()n+1 

9 


Av-^(Hn+1un+1) 


n 


(Hv)n+1  = {Hv)n  - 2At±—  Hnvn+1 

\J9o 

a21 

- 2At^—Hn+1vn+1  + 2At(G3  - G2)n 

\J9o 


+ 2AI-  B*  3 


iT Si)  Tbr) ) 


(Hn)2  5(7 


n+1 


Av-^{Hn+1vn+1) 


(2.106) 


(2.107) 
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2.9.5  Salinity  Scheme  (Two-Time- Level  Scheme) 

The  finite  difference  equation  for  salinity  are  similar  to  the  internal  mode  solution. 
The  vertical  diffusion  terms  are  treated  implicitly,  while  the  horizontal  diffusion  and 
advection  terms  are  treated  explicitly.  The  salinity  equation  in  the  (£,  r,,  a)  coordinate 
system  is  shown  in  the  following  (Sheng,  1986) 


8HS  B.8(  8S} 


dt 


HSCv  da  \ v da 


(2.108) 


where  G is  advection  plus  horizontal  diffusion  terms,  and 


A 

c _ VR 

^ Cv 


D 


VR 


C _ AHr 

bcH  ~ D 


G = 


Hr 

Ro 

\fdo 


+ |-( HvS + v&^(Hu,S) 


+ ^(g»™+29»d2S 


+ 9 


22 


d2S ' 


SCH  V * didr,  ' * dr,2 ) 

Using  i , j,  and  k as  grid  indices  in  the  £-,  //-,  and  cr-directions,  respectively,  the 
finite  difference  equation  for  salinity  can  be  written  as  (Sheng  and  Choi,  1989) 


H'+'SX,  - H'Sf/j,  = 


At  R. 


a 


+ ^(HvSJgiy  + 


(2.109) 


+ 


Ev(At)  1 
Hn+1SCv  Aak 


D 


A aA 


{sut  - sr') 


D,- 

A<7_ 


(sr1  - ss) 


+ 


HE„-At(  nd2S  , „ 12d2S 

a 9 + 2^f 


'CH 


di2 


didr, 


22  d2Ss 


dr,2 


Dividing  by  Hn+1  gives 


Ev  ■ At 


{Hn+1)2SCv  ■ A ak 


(Dy-  on+l  I Dy+  rrn+l  ) 

A<r_  •J.*-1  + a a+i’j’k+1) 


(2.110) 
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+ 


1 + 


Ev  ■ At  ( Dv+  Dv- 
(Hn+1)2SCvAak  \A^+A^: 


Cn+ 1 


II"  Eh  ■ At  ( ufj‘S 

H”+>  sc„  ^ ae 


H ' I!1"  W +2<712  ~ ^ ■ ~22 


^ 2od2S\" 

+ <72 


9^7/ 


dr/2 


Rearranging  Eq.  (2.110)  gives 


+ [l+T/iK-l  + ^+iiJSS-l^wiSaw.  (2-111) 


Tin 

12  on  i nn 

- Hn+1  z i,j,k  + '-'k 


This  is  applied  for  1 < k < km:  where 


7 k = 

EvAt 

(2.112) 

(H";')2Sc„  Aak 

ak+ 1/2  = 

Dv+ 
A cr+ 

(2.113) 

ak- 1/2  = 

i 1 
^ < 

(2.114) 

G*  = 

advection  + diffusion 

(2.115) 

For  k = 1,  i.e.,  the  lowest  a- level: 


1 + 7i«3/2]  - Jia3/2S?+} 


Hn 

Hn+l 


nn 


+ G? 


(2.116) 
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For  k = km,  i.e.,  the  upper- most  cr-level: 


+ [l  + lkm-l 


akm—  j 


r>n+X 


(2.117) 


Hn 


an  I rm 

^ i.i.lcm  ■ ^ leu 


J]n+ 1 i,j,km  ' ^ km 

Equations  (2.111),  (2.116),  (2.117),  can  be  written  as  the  following  tridiagonal  system 
of  algebraic  equations 

A1(K)  * S(K  - 1)  + A2{K)  * S(K)  + A3(K)  * S{K  + 1)  = A4(K)  (2.118) 


where  ,41,  ,42,  ,43,  and  ,44  are  defined  as 
A\{K) 


- - 


0 K = 1 

7*a*-i/2  1 < K < KM 


A2(K)  = 


A3(K)  = 


1 + 7ia3/2  K — 1 

1 + 1k(oikm-l/2  + (*k+ 1/2)  1 < K < KM 

1 T 7fcm— l/l^km  — 1/2  K KM 

~lk(*k+i/2  1 < I<  < KM 
0 A"  = KM 


(2.119) 


/U(A-)  = J^S^  + G?  1 < K < KM 


The  above  tridiagonal  system  of  equations  given  by  Equation  (2.118)  is  efficiently 
solved  by  a tridiagonal  solver  first  developed  by  Thomas  (1954),  which  is  often  called 
“The  Thomas  Algorithm”. 

2.9.6  Advection  Terms 

1.  Central  Difference  Scheme 

The  central  difference  option  has  the  advantage  of  second-order  accuracy,  but  a 
disadvantage  of  numerical  dispersion.  When  a large  concentration  gradient  exists  in 
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the  model  domain  with  insufficient  grid  points  across  the  gradient,  negative  concen- 


trations may  sometimes  develop  in  the  simulation.  The  finite- difference  form  is 


_1 

M [ 


{.HSy/gZ)ite  + (HSy/gdi+ijt 


Ui+l,j,k 


{HSJgdiM  + (HSy/gZi- itj,k 


ui,j,k 


(2.120) 


d_ 

dr] 


{HvSyfa)  = 


(HSy/g^v)+  - ( HSy/g~0v )_ 
At] 

1 

At] 


(HS^)i}jtk  + {HSy/gl) 

«,j+u 


' vi,j+l,k 


{H S y/9o)i,j,k  ~f  {H Sy/9o)i,j— l,k 


• ViJ,k 


(2.121) 


d^(Hu)S)  — ^ Uk  • ^(Si  J,k  + Si,j,k+ 1)  Uk- 1 • 2 {S*j,k  + Sij,k-1  )|  (2.122) 

2.  Upwind  Difference  Scheme 

The  upwind  difference  option  has  the  advantage  of  increased  stability  versus  the 
central  difference  option,  but  at  a cost  of  increased  numerical  diffusion.  The  upwind 
difference  option  guarantees  positive  salinity  concentrations.  The  finite-difference 
form  is 


= (S+uf+,J<k  - S.uJjj,)/ A£ 


(2.123) 


where  uT  = Hu^/g^. 

if  uI+i,j,k  > 0 


S+  = S, 


*,i,k 


(2.124) 


if  uI+i,j,k  < 0 


(2.125) 


Thus 
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S+uI+i,j,k  — 2 (uT+i,j,k  + lu  * + l)i> ^l)  Si,j,k  (2.126) 

"h  2 (u*+i ij’ifc  — lu«+i,j,fcl) 

Likewise, 

S- uTjj,  = 1 («^  + |«&i|)  Si-u*  (2.127) 

+ \ (uli,k  - i“Lii) 

3.  Combined  Central- Upwind  Difference  Option 

The  combined  central  and  upwind  option  gives  the  advantages  of  central  and 
upwind  options  together,  but  minimizes  the  disadvantages.  The  combined  option  is, 
however,  more  computationally  demanding.  The  finite-difference  form  is 

y((HuS^)  = (S+uf+w,t  - S.uJjj,)/ A(  (2.128) 

where  uT  = Hu^gl. 


if  uL,k  > 0 and  Sijjt  > Si+ij,k  : 

C T Ci+lJ,k 

S+~  2 

(2.129) 

if  ujj<k  > 0 and  SiJ>k  < Si+1<j<k  : 

S+  — Sij,k 

(2.130) 

if  ujhk  < 0 and  StJ,k  > Si+1jtk  : 

$+  = Si+l,j,k 

(2.131) 

d Ui,j,k  < d and  Sijje  < Si- klj,k  : 

C Sijje  + Si+itjtk 

S+~  2 

(2.132) 

4.  Higher  Order  Advective  Schemes 

To  improve  the  accuracy  of  advective  terms  in  the  numerical  model,  it  is  possible 
to  use  higher  order  advective  schemes  which  have  formal  truncation  errors  that  are 
third  order  or  fourth  order  in  Ax.  For  example,  Sheng  (1983)  applied  the  Flux- 
Corrected-Transport  (FCT)  scheme  (Zalesak,  1979)  to  simulate  the  advection  and 
diffusion  of  salinity  and  suspended  sediments  by  wind-driven  currents,  and  compared 
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the  performance  of  the  FCT  scheme  versus  the  three  advective  schemes  discussed 
above.  The  performance  of  the  combined  central  upwind  option  was  found  to  be 
quite  comparable  to  that  of  the  FCT  scheme.  More  recently,  Johnson  et  al.  (1989) 
applied  the  QUICKEST  scheme  developed  by  Leonard  (1979)  to  the  curvilinear  grid, 
and  obtained  improved  results  when  the  2-grid  version  of  CH3D  was  used.  Sheng  et 
al.  (1989)  found  that,  when  a cr-grid  model  is  used,  higher-order  advective  schemes 
should  be  avoided  in  regions  of  sharp  bathymetric  gradient  when  insufficient  grid 
points  exist  to  resolve  the  sharp  gradient. 

The  higher-order  schemes  have  been  primarily  applied  to  solve  the  advection 
and  diffusion  equations  but  not  to  solve  the  Navier-Stokes  equations.  Several  other 
higher-order  schemes  can  be  used  to  improve  the  accuracy  of  advective  terms.  The 
application  of  higher-order  advective  schemes  requires  additional  programming  and 
computational  effort,  due  to  the  fact  that  a differencing  formula  involving  four  or  five 
points  is  needed  to  represent  the  first-order  derivative  term. 

The  QUICKEST  scheme  used  for  this  study  has  third  order  accuracy.  This 
method  uses  cubic  upstream- weighted  interpolation  with  five  points  of  i — 2,  i — 1 , i,  i + 
l,and  i -f  2.  Because  this  is  a five-point  interpolation  the  fluxes  near  boundaries  must 
be  calculated  using  lower  order  schemes. 

The  finite  difference  equations  for  advection  terms  using  the  QUICKEST  scheme 
are  in  £- direction, 


d_ 

d( 


(HuSy/g~0)  = 


+ 


|uH-l,j,*:|)  2 + Si+l,j,k) 

— (l  — («t+lj,fcAQ  j (*Si+2)j)k  — 2Si+\jtk  + Si,j,k) 
2Ui+l,j,kAt  (Si+ ij}k  ~ Sitjtk) 

(ui+l,j,k  + lui+l,j,fc|)  2 + ‘S|'+l,j',fc) 

— (l  — (ui+i,j,fcAt)  ^ (<S't+ij,Af  ~ 2 Sij'k  + 
2Ui+l,j,k&t  (*^i+l ,j,k  ~ $ i,j,k ) 2 y/9°i+i ,j,u 
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(ui,j,k  |ut'j,fc|)  2 (Si- lj.fc  + Si,j, k)  (2.133) 

- g (!  ~ (ui,j,k&t)2)  ($i+i,j,k  - 2 Sijtfc  + Si-ij,*) 

- 2Ui,i,k^  (Si,j,k  - Si- U,fc) 

+ + |u«,j,*|)  ^ (Si—l,j,k  + -S'i.j.fc) 

- g (l  - (u.\j,*Af)2)  (5,-,^  - 25,_i,j,fc  + 5,_2,j,fc) 

- 9 u«j,fcA<  (Sij'k  — Si-i'j'k)  — 

z z y/9°i,jxni,j£ 


and  in  the  //-direction, 


^ (HvSy/g^) 


(vi+i,j,k  lu»+i,j,fc|)  2 ( Si,j,k  + Si+ lj.fc) 

g (l  - (u.+l,j,*A<)2)  (5't+2,j,fc  - 2S't+iJ)fc  + Sij'k) 

2U«+i,i,fcAt  (S't+iIj,fc  - 

(U*+l,j,fc  d"  lui+lj,fc|)  2 (Si,j,k  + ‘S't+l.j,*) 

— (l  — (n,+x,j,fcA<)  ) (S',+i)j)fc  — ZSij'k  + Si-ij'k) 

2 v«+i,j,fcAt  (Si+ij'k  — Sitjfk)  -y/g°t+l  j V H{+ijy 

(vi,j,k  ~ luij,*|)  2 (Si-1, j,k  + 5'*,i,fc ) (2.134) 

- (l  — (utj,fcAi)  ) (•S'i+ij^  — 25'ij)*:  -f  •S't-ij,*;) 

-u,ijifcAt  (S.-j,*  - Si-u,*) 

(vi,j,k  + lU«,j,fc|)  2 (Si  — l,j,k  + Sij'k) 
g (1  ~ (u.\j,fcA*)2)  (5,-j,*  - 2Si-itjtk  + Si-2, j,k) 

nVi,j,k&t  (Sij'k  ~ Si-lJtk)  9 y/QoiJ.V  Hi,J,V  ~~P- rr 

Z Z y/9°i,)Xni,jA 


2.9.7  Horizontal  Diffusion  Terms 


The  finite- difference  form  of  the  horizontal  gradient  in  the  horizontal  diffusion 
terms  is  the  standard  second-order  divided  difference.  There  is  no  need  to  treat  the 
horizontal  diffusion  terms  implicitly  because  the  time  step  constraint  associated  with 
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them  are  usually  much  larger  than  that  allowed  by  the  advection  terms. 


d2S  n Si+itjtk  — 2 Sijtk  + Si-ij'k 


ae 


= 9 


(AO2 


(2.135) 

(2.136) 


,12 


- 9u  [(‘S'.+i/2,i+i/2,fc  - Si+i/2,j-i/2,k)  (2.137) 
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2.9.8  Numerical  Stability 


(2.138) 


Since  the  numerical  scheme  used  here  is  not  fully  implicit,  the  time  step  of  numer- 
ical computation  cannot  exceed  the  limits  associated  with  the  advection  terms,  the 
Coriolis  terms,  the  baroclinic  pressure  gradient  terms,  and  the  horizontal  diffusion 
terms 


( At)advection 
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z> 

(At)  Coriolis 
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\[^p  9 Umax  ( aJ  + Ay  )m“* 

(2.139) 
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H 


1 I 1 

(Ax^fA^p- 


If  the  allowable  internal  time  step  according  to  the  above  equations  is  much  (10 
to  100  times)  larger  than  the  At  allowed  by  the  surface  gravity  wave  as  shown  in 
Equation  (2.139),  then  it  is  feasible  to  use  a small  time  step  for  the  internal  mode 
(A ti).  If  the  difference  between  the  allowable  internal  time  step  is  only  a few  times 
larger  than  the  external  step,  then  it  is  better  to  simply  use  the  same  A te  and  At,. 
The  selection  of  time  step  is  still  very  much  an  “art”  because  all  the  time  step  limits 
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are  based  on  stability  analysis  which  is  strictly  valid  only  in  a linear  system  without 
considering  the  effects  of  nonlinearity  and  boundary  conditions. 

2.10  Turbulence  Model  for  Mass  and  Momentum  Exchange  Coefficients 

Various  models  are  available  to  parameterize  the  vertical  turbulent  mixing  of 
momentum,  heat,  and  salinity  in  the  water  column.  For  reviews  of  turbulence  models, 
the  readers  are  referred  to  Lewellen  (1977),  Rodi  (1980),  Mellor  and  Yamada(1982), 
and  Sheng  (1986).  The  starting  point  of  a turbulence  model  is  to  derive  the  mean  flow 
equations  from  the  conservation  equations  for  mass,  momentum,  energy,  and  species. 
For  incompressible  flows,  conservation  equations  can  be  written  in  tensor  notation  as 


7T  = ° 

OX; 


(2.140) 


du{ 


dui  _ 1 dp  (, p - po) 


dt  3 dxj  po  dxi  p0 


d ( dui 


‘2‘£ijk^Jjuk  To  [ v 


dxj  y dxj 


d<t>  , d<t>  , d f d<t> 
H7  + uJinr  = +—'k- 


(2.141) 


(2.142) 


dt  3 dxj  dxj  \ dxj) 
where  (j)  is  either  density  or  temperature  or  salinity,  p is  density,  po  is  reference  density, 
£{jk  is  the  permutation  tensor,  f lj  is  the  rotational  speed  of  earth,  and  v and  k are 
the  molecular  viscosity  and  diffusivity. 

When  the  variables  are  separated  into  mean  and  fluctuating  quantities,  the  mean 
flow  equations  can  be  written  as 


(2.143) 


dui  , duj  = du\u'3  1 dp  (p  - Po) 

dt  3 dxj  dxj  p0  dxi  pQ 

o o — 9 ( dui\ 

2cijt(ljUk  + g—  (^j/— J 

d}  _ du'rf  d ( t d4>\ 

dt  3dXj  dXj  + dXj  [ dxj 


(2.144) 


(2.145) 
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where  u-,  uj,  and  are  the  mean  velocity  components,  u\  and  u'  are  the  fluctuating 
velocity  components,  <j>  and  <j>'  are  the  mean  and  fluctuating  temperature  or  salinity, 
and  p is  the  mean  pressure.  This  system  of  equations  is  not  closed  because  of  the 
Reynolds  stress  terms. 

Basically,  turbulence  models  can  be  classified  as  eddy  viscosity  models  and  second- 
order  closure  models.  The  eddy  viscosity  models  assume  the  turbulent  Reynolds 
stresses  to  be  the  products  of  mean  velocity  gradients  and  “eddy  viscosities,”  and 
the  turbulent  heat  and  mass  fluxes  to  be  the  products  of  mean  temperature  and 
concentration  gradients  and  “eddy  diffusitivities” 


— , (dui  duj\  q\ 

(2.146) 

OXi 

(2.147) 

where  q 2 = u'u'-,  and  At  and  Kt  are  turbulent  eddy  viscosity  and  diffusivity,  respec- 
tively. 

The  eddy  coefficients  are  assumed  to  be  certain  known  functions  of,  for  example, 
mean  flow  parameters  (Munk  and  Anderson,  1948),  depths  and  wind  speeds.  These 
functions,  however,  are  difficult  to  determine  and  are  often  ad-hoc  in  nature,  and 
hence  may  have  to  be  adjusted  for  application  to  a new  site  or  a new  flow  situation. 
Obviously,  this  limits  the  “predictability”  of  the  hydrodynamics  model. 

The  second-order  closure  models,  on  the  other  hand,  resolve  the  dynamics  of 
turbulence  by  including  the  differential  transport  equations  for  the  turbulence  vari- 
ables; i.e.,  the  second-order  correlations  (— and  W).  In  the  most  compli- 
cated case,  i.e.,  the  Reynolds  stress  model,  differential  transport  equations  are  solved 
(Lewellen,  1977;  Sheng,1982;  Sheng  and  Villaret,  1989).  The  time  averaged  mo- 
mentum equations  are  subtracted  from  the  instantaneous  Navier-Stokes  equations  to 
produce  the  primed-variable  equations.  The  resulting  equation  for  the  i-component 
is  multiplied  by  u'-,  and  j-  component  equation  by  u\.  Adding  two  equations  and 
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averaging  by  time,  gives  the  ifjifj— equation. 


duty 

dt 
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+ Uj  dxj  - ^ dx3 
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dxk 


dx2 
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dx'  dxj 

where  x,  are  coordinate  axes,  t is  time,  (iq,  uj,  u^)  are  the  mean  velocity  components, 
(uj-,«j,  ujb)  are  the  fluctuating  velocity  components,  (j)  and  <f>'  are  the  mean  and  fluc- 
tuating values  of  density  or  temperature  or  salinity,  eijk  is  the  alternating  tensor,  and 
fl  is  earth  rotation. 

The  above  system  of  equations  is  not  closed  since  each  equation  contains  unknown 
terms.  Following  the  second-order  closure  model  described  in  Lewellen  (1977)  and 
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duW.  du'V 

= 


dt 


Sheng  (1982),  the  above  equations  are  replaced  by  the  following  equations. 

(Evolution) 
(Production) 

(Buoyancy) 
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(2.153) 


where  Sy  is  the  Kronecker  delta,  q = (uju-)1/2  is  the  total  rms  fluctuating  velocity, 
and  A is  the  turbulence  macroscale  which  is  a measure  of  the  average  turbulent  eddy 
size.  A total  of  five  model  coefficients:  a,  A,  5,  wc,  and  s are  contained  in  the  above 
equations.  By  matching  model  results  to  data  from  critical  laboratory  experiments 
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where  only  one  or  two  of  the  turbulent  transport  processes  are  important,  these  model 
constants  were  found  to  be  (Lewellen,  1977)  a = 2.5,  A = 0.75,6  = 0.125,  vc  = 0.3 
and  s = 1.8.  An  additional  equation  for  A is  needed  to  close  the  system. 

The  present  study  tested  several  methods  for  calculating  the  vertical  turbulent 
eddy  coefficients:  (1)  constant  eddy  coefficients,  (2)  Munk- Anderson  type  eddy  coef- 
ficients, (3)  a simplified  second-order  closure  model  of  turbulent  transport,  i.e.,  the 
equilibrium  closure  model,  and  (4)  another  simplified  second-order  closure  model,  i.e., 
the  turbulent  kinetic  energy  (TKE)  closure  model.  These  methods  are  discussed  in 
the  following. 

• Constant  Eddy  Coefficients 


Constant  eddy  coefficients  are  generally  used  in  preliminary  model  runs.  The 
rule  of  thumb  is  that  the  vertical  eddy  coefficients  should  increase  with  increasing 
wind  stress  and/or  tidal  amplitude.  Coefficients  can  be  adjusted  until  a reasonable  fit 
exists  between  model  results  and  data.  In  the  absence  of  data,  however,  it  is  difficult 
to  adjust  the  constant  eddy  coefficients  to  give  good  results. 

• Munk-Anderson  Type  Eddy  Coefficients 

In  non-stratified  flow  situations,  a variable  eddy  coefficient  of  the  following  form 
is  used 


A 


Vo 


ZrUrAl  17 M 2 (dv\ 

H [da)  +{da) 


(2.154) 


where  u and  v are  dimensionless  horizontal  velocities,  the  turbulence  macroscale  A0 
is  assumed  to  be  a linear  function  of  a , increasing  with  distance  above  the  bottom  or 
below  the  free  surface,  and  with  its  peak  value  at  mid-depth  not  exceeding  a certain 
fraction  of  the  local  depth.  In  the  presence  of  strong  waves,  turbulent  mixing  in 
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the  upper  layers  may  be  significantly  enhanced.  In  such  a case,  the  length  scale  A0 
throughout  the  upper  layers  may  be  assumed  to  be  equal  to  the  maximum  value  at 
mid-depth. 

In  stratified  flow  situations,  the  effect  of  buoyancy  on  eddy  coefficients  can  be 
parameterized  in  terms  of  several  Richardson  number-dependent  stability  functions 

Av  = Avo<f>\(Ri)',  Kv  = Kvo<f>2(Ri);  Dv  = DV0(j)2{Ri)  (2.155) 

where 

_ —gHZre  dp 
Kl~  U?(l+ep)da 

where  Avo,  Kvo,  and  Dvo  are  eddy  coefficients  in  the  absence  of  any  density  stratifica- 
tion, and  (f>i  and  <j)2  are  stability  functions  traditionally  assumed  to  be  of  the  following 
forms 

<h  = {l  + a1Jtt)A;  fa  = (l  +a2Ri)lh  (2.157) 

As  discussed  by  Sheng  (1983),  there  exist  great  discrepancies  among  the  many 
existing  empirical  forms  of  the  stability  functions.  To  determine  the  coefficients  in 
these  stability  functions,  one  needs  sufficient  data  to  achieve  the  “best  fit”  between 
modeled  and  measured  results.  As  such,  the  specific  formulae  and  coefficients  will 
vary  with  the  problem,  the  environment,  and  the  available  data.  Munk  and  Anderson 
(1948)  used  the  following  formula  in  a study  of  the  marine  thermocline 

<h  = (1  + 10  Ri)~1/2-,  <t>2  = (1  + 3.33  Ri)~3/2  (2.158) 

Others  (Kent  and  Pritchard,  1959;  Blumberg,  1975;  Bowden  and  Hamilton,  1975) 
used  <j) i and  <f>2  in  the  form  of  Equation  (2.157)  but  with  coefficients  different  from 
those  in  Equation  (2.158).  Not  only  the  form  of  the  stability  function  may  vary 
with  the  problem,  but  different  Richardson  numbers  may  also  be  used  for  different 
types  of  problems.  For  example,  the  formation  and  deepening  of  the  thermocline  in 


'duV  /dv\2 
,daj  + [ da j 


(2.156) 
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a relatively  shallow  basin  depends  strongly  on  the  relative  importance  of  wind  stress 
and  heat  flux  at  the  free  surface.  In  such  a case,  the  following  Richardson  number 
could  be  used 


n2gHZra2e  dp 

r\  7 — 

U2ul{  1 + ep)  da 

where  ac  is  the  von-Karman  constant,  and  u*  is  the  dimensionless  friction  velocity  at 
the  free  surface. 

• A Simplified  Second-Order  Closure  Model:  The  Equilibrium  Closure  Model 

An  equilibrium  closure  model  (Sheng,  1985;  Sheng  and  Chiu,  1986;  Sheng  et  al., 
1989)  has  been  used  to  compute  the  vertical  turbulence.  As  discussed  before,  the  local 
equilibrium  condition  is  valid  when  the  time  scale  of  mean  flow  is  much  larger  than 
that  for  turbulence  and  when  turbulence  varies  little  over  the  turbulence  macroscale. 
In  this  case,  the  evolution  terms  and  the  diffusion  terms  in  Equations  (2.151  )-(2. 153) 
become  negligible  with  respect  to  other  terms.  The  equilibrium  closure  model  is 
significantly  simpler  than  the  complete  second-order  closure  model  (Reynolds  stress 
model)  and  has  been  found  to  give  very  good  results  in  mean  flow,  salinity  and 
temperature  with  little  or  no  tuning  of  model  coefficients. 

In  addition  to  the  mean  flow  equations,  a set  of  algebraic  equations  are  solved  for 
the  second-order  correlation  quantities  to  obtain  the  stability  functions  <f>\  and  <^2  in 
terms  of  the  mean  flow  variables.  These  equations,  when  written  in  dimensional  and 
tensor  forms,  are 
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where  the  subscripts  i,j,k  can  take  on  the  values  of  1,  2,  or  3.  Hence  Equation 
(2.160)  represents  six  equations  and  Equation  (2.161)  represents  three  equations  in 
general.  A total  of  five  model  coefficients  are  contained  in  the  above  equations. 
These  coefficients  were  determined  from  comparing  model  results  with  data  from 
critical  laboratory  experiments  where  only  one  or  two  of  the  turbulent  transport 
processes  are  dominant.  These  coefficients  remained  “invariant”  in  application  of 
the  equilibrium  closure  model,  the  kinetic  energy  closure  model  (Sheng  and  Villaret, 
1989)  and  the  Reynolds  stress  model  (Sheng,  1982;  Sheng,  1984;  Sheng  and  Villaret, 
1989).  A graphical  comparison  of  this  formulation  versus  some  empirical  formulae  is 
shown  in  Figure  2.10. 

As  shown  in  Sheng  et  al.  (1989),  q 2 can  be  determined  from  the  following  dimen- 
sional equations  when  mean  flow  variables  are  known 


3 A2b2sQ4  + A[(bs  + 3b  + 7b2s)Ri-  Abs(l  -2b)]Q2  (2.163) 

+ b(s  + 3 + 4bs)Ri2  + (bs-A)(l  -2b)Ri  = 0 

where  A,  b , and  s are  model  constants,  Ri  is  the  Richardson  number  and 


q = Q a 


\ 


(du\‘  fdv\2 

(si)  + y 


(2.164) 


Once  q2  are  computed,  the  vertical  eddy  coefficients  can  be  computed  from 


<t>2(FH) 
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(a) 


(b) 


Figure  2.6:  Comparison  of  stability  functions  for  (a)  some  epirical 
formulae  versus  (b)  the  simplified  second-order  closure  model  (From 
Sheng,  1983). 
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Av 

Kv 

where 

w = Ri/(AQ 2)  (2.167) 

w = w/(l—w/bs)  (2.168) 

W = 3(I^)?2  (2'169) 

In  the  complete  second-order  closure  model,  a differential  transport  equation  for 
the  turbulent  macroscale  A is  derived.  The  equation  is 


A + ww'w' 
A — wq2 
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However,  the  A equation  contains  four  model  coefficients  that  must  be  calibrated 
with  experimental  data.  For  ease  of  application  but  without  loss  of  generality,  the 
turbulent  macroscale  A is  assumed  to  satisfy  a number  of  integral  constraints.  First 
of  all,  A is  assumed  to  be  a linear  function  of  vertical  distance  immediately  above 
the  bottom  or  below  the  free  surface.  In  addition,  the  turbulent  macroscale  A must 
satisfy  the  following  relationships  (Sheng  and  Chiu,  1986;  Sheng  et  al. , 1989) 


dA 

dz 


< 0.65 


(2.171) 


A < Ci-H 


(2.172) 
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(2.173) 


A < Cx-Hr 


A < C2-Sq  2 


(2.174) 


(2.175) 


where  Cx  is  a number  usually  within  the  range  of  0.1  to  0.25;  H is  the  total  depth;  Hp  is 
the  depth  of  the  pycnocline;  C2,  which  is  between  0.1  and  0.25,  is  the  fractional  cut-off 


in  Lake  Okeechobee  (Sheng  et  al.,  1989b)  gave  good  results  with  little  or  no  tuning 
of  model  coefficients. 

• Turbulent  Kinetic  Energy  (TKE)  Closure  Model 

A simplified  TKE  closure  model  (Sheng  and  Villaret,  1989)  has  been  used  to  compute 
the  vertical  turbulence.  When  the  time  scale  of  the  mean  flow  evolution  is  much 
greater  than  the  time  scale  of  turbulence  (A /q),  the  evolution  terms  and  the  diffusion 


limitation  of  turbulent  macroscale  based  on  8q 2,  the  spread  of  turbulence  determined 
from  the  turbulent  kinetic  energy  ( q 2)  profile;  and  N is  the  Brunt-Vaisala  frequency 
defined  as 


(2.176) 


Although  the  simplified  second-order  closure  model  presented  above  is  strictly 
valid  when  the  turbulence  time  scale  (A /q)  is  much  less  than  the  mean  flow  time 
scale  and  when  turbulence  does  not  change  rapidly  over  A,  it  has  been  found  to  be 
quite  successful  in  simulating  vertical  flow  structures  in  estuarine  and  coastal  waters. 
Recent  simulations  of  wind-induced  mixing  of  salinity  in  Chesapeake  Bay  (Sheng  et 
al-,  1989a;  Johnson  et  al. , 1989)  and  wind-driven  circulation  and  sediment  transport 
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terms  in  the  second-order  closure  model  become  negligible  with  respect  to  other  terms. 
Instead  of  calculating  q 2 by  algebraic  equations,  q 2 is  obtained  in  a dynamic  equation 
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In  the  one-dimensional  vertical  case  of  the  <x-grid,  the  equation  is 
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(2.178) 


Applying  the  rules  in  section  2.3,  nondimensionalization  of  the  above  equation  gives 
(all  the  variables  are  now  nondimensional) 


dq2  Ro  dHuq 2 Rovdojq 2 1 d . dq2 . 

nEvAV/du^2  , 2EVKV  dp  n q 3 
H*  [da’  + Fr2DH  da  v 4 A 


(2.179) 


where  and  Tf,,  are  nondimensionalized  with  Aur,  /?o„  = f/r//Zr  is  a vertical  type 
Rossby  number,  and  Ev  = Avr/fZ2. 

Again,  integral  constraints  for  A are  used  as  the  equilibrium  closure  model. 


CHAPTER  3 

MODEL  ANALYTICAL  TEST:  COMPARISON  BETWEEN  MODEL  RESULTS 

AND  ANALYTIC  SOLUTIONS 

When  numerical  models  are  used  to  simulate  circulation  and  transport  in  realistic 
estuarine  environments,  it  is  common  to  find  discrepancies  between  numerical  results 
and  data.  The  discrepancies  often  result  from  a variety  of  errors:  errors  associated 
with  measurement,  errors  associated  with  the  formulation  of  the  physical  problem 
(e.g.,  dimensionality,  turbulence  closure  scheme,  moving  boundary  scheme,  etc.),  er- 
rors associated  with  the  formulation  of  the  numerical  problem  (e.g.,  finite-difference 
schemes  or  finite-element  schemes  for  various  terms  of  the  governing  equations),  and 
errors  associated  with  the  computer  and  programming  (e.g.,  computer  round-off  error 
and  programming  error).  Before  the  model  is  applied  to  a realistic  estuary,  it  is  nec- 
essary to  test  and  improve  the  numerical  accuracy  of  the  various  numerical  schemes 
used  in  the  model.  This  chapter  presents  a study  on  the  numerical  accuracy  of  the 
three-dimensional  hydrodynamic  model  by  comparing  model  results  with  a number  of 
analytic  solutions  describing  the  circulation  in  idealized  basins  forced  by  a constant 
surface  slope,  seiche,  tide,  wind,  and  density  gradients. 

3.1  Forced  Constant  Surface  Slope 

The  first  test  considers  the  steady-state  response  of  a water  body  in  a rectangular 
basin  forced  by  constant  surface  slope  with  uniform  density.  At  steady  state,  the 
surface  slope  balances  the  bottom  stress.  Assuming  one-dimensional  flow  with  the 
quadratic  bottom  friction,  the  governing  equation  is  given  as 
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~9hiL = Cd  i u i u 


(3.1) 


with  boundary  conditions  ( = a at  i = 0,  and  £ = -a  at  x = l (Figure  3.1), 
where  u is  the  vertically-averaged  velocity,  £ is  surface  elevation  above  mean  sea  level, 
h is  water  depth,  g is  gravitational  acceleration  (=  981  cm/s2),  l is  the  basin  length, 
and  Cd  is  the  bottom  friction  factor.  A non-uniform  rectangular  grid  as  shown  in 
Figure  3.1  was  used.  The  basin  is  open  at  the  left  and  right  ends.  Model  simulation 
was  started  with  zero  elevation  and  velocity  initially.  Model  tests  were  performed 
for  the  cases  of  linear  bottom  friction  and  quadratic  bottom  friction.  The  following 
parameter  values  were  used 


a = 25  cm 


/ = basin  length 
basin  width 
h 

Cd 


= 62  km 
= 14  km 
= 10  m 


(3.2) 


0.003  for  quadratic  bottom  friction 
0.3  for  linear  bottom  friction 


According  to  Equation  (3.1),  the  steady  state  velocity  was  found  to  be  u = 
52.17  cm /sec  for  the  quadratic  bottom  friction  case,  and  u = 27.2  cm/ sec  for  the 
linear  bottom  friction  case. 

Using  a time  step  of  360  seconds,  model  results  reached  steady  state  in  about  300 
time  steps.  The  steady-state  vertically-averaged  velocities  in  the  x-y  plane  as  shown 
in  Figure  3.2  agree  well  with  the  analytical  results.  Figure  3.3  shows  the  velocity 
distribution  in  a basin  which  is  rotated  in  the  x-y  plane  by  30°  in  the  counterclockwise 
direction.  The  results  appear  to  be  the  same.  The  model  was  then  tested  using  an 
irregular  curvilinear  grid  shown  in  Figure  3.4.  The  same  steady  velocity  is  reproduced 
(Figure  3.5). 
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Figure  3.1:  A non-uniform  rectangular  grid  system  for  testing  the  forced  surface  slope 
test. 


50  cm/sec 


Figure  3.2:  Simulated  steady-state  vertically-averaged  velocity  field  in  a rectangular 
basin  forced  by  a constant  surface  slope.  Horizontal  grid  is  shown  in  Figure  3.1. 
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50  cm/sec 


Figure  3.3:  Same  as  Figure  3.2,  except  that  the  basin  is  rotated  by  30°  in  the  coun- 
terclockwise direction. 
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Figure  3.4:  An  irregular  grid  system  for  rectangular  basin  with  irregular  grids 
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50  cm/sec 


Figure  3.5:  Same  as  Figure  3.2,  except  that  the  irregular  curvilinear  grid  as  shown  in 
Figure  3.4  is  used. 
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3.2  Seiche  Test 

The  second  test  was  designed  to  test  the  model’s  ability  to  simulate  a seiche 
(standing  wave)  in  a closed  basin.  The  local  acceleration  term  balances  the  pressure 
term.  Neglecting  vertical  and  horizontal  diffusion,  bottom  and  side-wall  friction, 
convective  acceleration,  and  Coriolis  acceleration,  the  equations  of  motion  reduce  to 
the  following  one-dimensional  equations 

du  d( 

dt  9 dx 

du  _ d( 
hdi  = ~dt 

with  boundary  conditions  u = 0 at  x = 0,  /,  where  u is  velocity,  ( is  surface  elevation, 
h is  water  depth,  g is  gravitational  acceleration  (=  981  cm/s2),  and  l is  the  basin 
length.  The  analytical  solutions  of  the  lowest  mode  for  the  above  two  equations  are 

( = a cos  ut  cos  kx  (3.5) 

gak  . . , . 

u = smutsmkx  (3.6) 

where  a is  the  wave  amplitude,  u>  is  the  circular  frequency,  and  k is  the  wave  number. 
The  grid  system  used  here  is  the  same  as  that  shown  Figure  3.1.  The  following  values 
were  used  for  this  test  case 


a 

l = basin  length 
basin  width 


15  cm 

62  km  = A/2 
14  km 
5 m 


(3.7) 


h 
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BASIN  LENGTH  FROM  ONE  END  (xll) 


Figure  3.6:  Comparison  between  simulated  surface  elevation  (diamond  symbol)  and 
analytic  solution  (solid  lines)  for  seiche  oscillation  in  a closed  basin  (T  is  seiche  period). 
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SEICHE  TEST 


BASIN  LENGTH  FROM  ONE  END  (xll) 


Figure  3.7:  Comparison  between  simulated  vertically- averaged  velocities  (diamond 
symbol)  and  analytic  solution  (solid  line)  for  seiche  oscillation  in  a closed  basin  (T  is 
seiche  period). 
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t=0  20  cm/sec 


Figure  3.8:  Maximum  velocities  for  seiche  case 
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The  wave  characteristics  are  calculated  as 

jfe  = 2tt/A  = 5.0671  • 10_7/cm 

T = A/ \Jgh  = 4.92  hours  (3-8) 

u = 27t/T  = 3.5488  • 10-4/s 

The  above  values  and  a time  step  of  300  seconds  were  used  in  the  numerical  simulation, 
using  the  vertically- integrated  version  of  the  three-dimensional  model.  Figures  3.6 
and  3.7  show  the  comparison  between  model  results  and  analytic  solution  after  the 
model  has  been  run  for  ten  tidal  cycles.  Figure  3.6  shows  the  surface  elevations 
distribution  along  the  channel  at  increments  of  one  eighth  tidal  period.  Figure  3.7 
shows  the  velocity  distribution  at  corresonding  times.  The  model  predictions  match 
the  analytic  solution  very  well  with  a maximum  error  of  3 % in  amplitude  of  both 
the  elevation  and  the  velocity.  Figure  3.8  shows  the  instantaneous  velocity  fields  at 
two  instants  of  time  when  the  velocity  is  maximum.  The  maximum  velocity  occurs 
at  the  center  and  zero  velocity  occurs  at  both  closed  ends. 

3.3  Tidal  Forcing 

One  of  the  most  common  applications  of  an  estuarine  hydrodynamic  model  is 
tidal  simulation.  Therefore,  before  applying  the  model  to  a real  field,  the  model  result 
should  be  compared  with  the  analytic  solution  describing  the  response  of  a rectangu- 
lar estuary  to  forced  oscillation  at  the  open  boundary.  Lynch  and  Gray  (1978)  derived 
the  analytic  solutions  for  a tidally  forced  estuary  with  a flat  bottom  and  a sloping 
bottom.  Neglecting  nonlinear,  Coriolis,  horizontal  diffusion,  and  ^-direction  terms 
and  assuming  linear  bottom  friction,  one  obtains  the  following  vertically-averaged 
equations 


* + 
at 


dhu 


= 0 


dx 


(3.9) 


du 

~di 
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(3.10) 


d(  ku 

+ 9di  + T~° 

with  the  following  boundary  conditions:  at  the  closed  end,  x/,  d(,/dx  = 0;  and  at 
the  open  end,  x0,  C = acosu;£,  where  u is  the  velocity,  ( is  surface  elevation,  g is 
gravitational  acceleration,  k is  the  linear  bottom  friction  coefficient,  a is  the  tidal 
amplitude,  u>  is  the  tidal  frequency,  and  h is  the  depth.  Generally  h is  defined  as 
h{x)  = H0xn , while  h with  a linear  slope  is  given  as 


L H0 -Hi 

h = x = sx 


(3.11) 


where  Ho  is  the  depth  at  the  open  end  (x  = 0),  Hi  is  the  depth  at  the  closed  end 
(x  = /),  l is  the  basin  length,  and  s is  the  bottom  slope.  Flat  bottom  and  linear 
bottom  slope  cases  are  considered  for  model  tests. 

The  solution  for  a flat  bottom  (n  = 0)  is  (Lynch  and  Gray,  1978) 


COM) 

u(x, t ) 


Re[ae 


iu ;tCOS(/?(x  - X/)) 

COS  (/?/) 


ftc[  ~iua  c%wt  sin(/3Qc  - */)) 


PH  o 


cos 


m 


(3.12) 

(3.13) 


where  a is  the  tidal  amplitude  at  the  open  end,  and 

0 = -iu>k/h  0'S 
9H0 


(3.14) 


The  solution  for  a linear  sloping  bottom  is  (Lynch  and  Gray,  1978) 


C(M)  = Re[{AJ0(2pyfi)  + BY0(2pyft))eiut]  (3.15) 

u(x,t)  = Yle[-^=[(— AJi(2fiy/x)  — BKi(2/?v^))]^~e'w(]  (3.16) 

\JX  ps 

where  Jq,J\,Yq,  and  Yi  are  Bessel  functions,  P is  redefined  with  s instead  of  H0,  and 


aY\(2Py/xi) 

[Jo(2Py/xo)Yi(2Py/x~i)  - Y0(2py/x^)Jl(2Pi/xi)\ 

-aJi{2py/xj) 

[Jo(2Py/^)Y1(2p^rl)  - Y0(2py/^)J1(2p^i)} 


(3.17) 

(3.18) 
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The  following  parameter  values  were  used  in  the  test  run 


a = 1 00  cm 
T = 600  sec 
k = 1 cm/s 

l = 6200  m 
(H0  + Hi) /2  = 350  m 
At  = 5 s 


(3.19) 


For  a flat  bottom,  Figure  3.9  shows  comparison  between  model  results  and  the 
analytic  solution  for  elevation  and  velocity  along  the  basin  length.  Maximum  error 
is  less  than  2%  for  elevation,  and  3%  for  velocity.  For  a sloping  bottom,  Figure  3.10 
shows  simulated  and  analytical  results  as  a function  of  bottom  slope  for  maximum 
surface  elevation  at  the  closed  end  and  for  velocity  at  the  open  end.  Simulated  velocity 
and  elevation  both  agree  well  with  the  analytic  solution. 

Another  tidal  forcing  test  was  done  with  the  annular  section  configuration  shown 
in  Figure  3.11.  Here  h is  defined  as  h(x)  = H0rn , and  /3  is  defined  as  /?2  = (u>2  — 
iu>k/h)/gH0.  Flow  is  required  to  be  tangent  to  the  solid  boundaries  at  r = r\.  A 
tidal  forcing  function  is  specified  at  r = r2.  The  analytic  solution  for  a flat  bottom 
(n  = 0)  is  (Lynch  and  Gray,  1978) 

C(r,t)  = Re[(A70(/?r)  + BY0^r))eiut]  (3.20) 

b(M)  = Re[(-/U1(^r)-By,(/?r))^-e'"']  (3.21) 

where  J0,  J\,Y0,  and  Vi  are  Bessel  functions,  and 

aY\{pri) 

Wr2)Yx{^)  -Y0^r2)J^rx)] 

-aJi(pr  i) 


B 


(3.23) 


ELEVATION  (Ml 
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TIDAL  FORCING 


OISTRNCE  FROM  THE  CLOSED  END  (X/L) 


Figure  3.9:  Comparison  between  simulated  surface  elevations  and  velocities  (*’s)  and 
analytic  solutions  for  a tidally  forced  flat-bottom  basin  (Elevation  and  velocity  are  in 
phase). 
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TIDAL  FORCING  FOR  SLOPING  BOTTOM 


LOG 10 (BOTTOM  SLOPE) 


Figure  3.10:  Comparison  of  simulated  maximum  surface  elevations  at  the  closed  end 
and  maximum  velocities  at  the  open  end  with  the  analytic  solution  for  a tidally  forced 
sloping  bottom  basin. 
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Figure  3.11:  Grid  system  for  the  annular  section 
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Figure  3.12:  Comparison  between  simulated  surface  elevations  and  velocities  (*’s) 
and  analytic  solutions  for  a tidally  forced  flat-bottom  annular  section. 
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50.0  CM/SEC 


Figure  3.13:  Maximum  ebb  and  flood  velocities  in  a tidally  forced  flat-bottom  annular 
section  at  the  end  of  tidal  10  cycles  run. 
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A sinusoidal  tidal  forcing  was  prescribed  along  the  outer  radius  for  the  model 
simulation.  The  following  values  were  used  in  the  test  run 

a — 100  cm 
T = 2.5  hours 

k = 1 cm/s  (3.24) 

/ = 63  km 
h = 100  m 
At  = 100  s 


The  model  was  run  up  to  10  cycles,  at  which  steady  state  conditions  were  obtained. 
Figure  3.12  shows  the  resulting  surface  elevation  and  velocity  at  maximum  conditions. 
Model  results  appear  to  agree  well  with  the  analytic  solution.  Maximum  ebb  and  flood 
velocities  are  plotted  in  Figure  3.13.  The  radial  symmetry  is  also  precisely  reproduced 
in  the  model  results. 

3.4  Wind  Forcing 


Three  simulations  of  wind  forcing  have  been  done.  The  first  test  was  conducted 
to  examine  the  effect  of  bottom  slope  on  surface  setup  in  a rectangular  basin  with  one 
open  end.  The  second  test  was  to  examine  the  vertical  profile  of  horizontal  velocity 
in  a closed  basin.  The  third  test  was  to  examine  the  surface  setup  in  bowtie-shaped 
estuary. 

The  analytic  solution  for  the  first  test  was  solved  by  Lynch  and  Gray  (1978)  by 
neglecting  the  convective  acceleration,  horizontal  friction,  Coriolis,  and  y-direction 
terms,  while  assuming  small  amplitude  surface  oscillation  and  using  a linearized  bot- 
tom friction.  Thus  the  vertically-averaged,  linearized  shallow  water  equations  are 
given  as 


dt 


+ 


dhu 


= 0 


dx 


(3.25) 
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du  d(  ku  tw 

dt  9 dx  h ph 


(3.26) 


with  the  boundary  condition  at  the  closed  end,  xj,  as 

_ Jw_ 
dx  pgh 


(3.27) 


and  at  the  open  end,  x0,  as 

C = 0 (3.28) 

where  u is  the  velocity,  ( is  surface  elevation,  g is  gravitational  acceleration,  k is  the 
linear  bottom  friction  coefficient,  tw  is  the  wind  stress,  p is  the  water  density,  and  h 
is  the  depth.  Generally  h is  defined  as  h(x)  = Haxn  while  h with  a linear  slope  is 
given  as 

h = — — -x  = sx  (3.29) 

where  Ho  is  the  depth  at  the  open  end,  Hi  is  the  depth  at  the  closed  end,  l is  the 
basin  length,  and  s is  the  bottom  slope. 

The  analytic  solution  for  the  steady  wind  setup  is 


C - rr  (x  xo)  for  n - 0 
pgtio 

(3.30) 

— _Jfl(l0g  x _ ]0g  x \ for  n = 1 

P9S 

(3.31) 

where  x0  is  the  coordinate  at  the  open  end,  while  x/  is  the  coordinate  at  the  closed 
end  (0  < X/  < x < xQ). 

For  the  second  test,  the  analytic  solution  for  the  vertical  velocity  profile  driven 
by  a constant  wind  stress  in  a closed  basin  can  be  readily  obtained.  Considering  only 
pressure  and  vertical  diffusion  terms,  the  governing  equations  are  given  as 


d(  Ad2u 
9 dx  ~Adz 2 

(3.32) 

[0 

/ udz  = 0 
J-h 

(3.33) 
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with  the  boundary  condition  at  the  surface  (z  = £)  being 


. du  Tw 

aTz  = J 


(3.34) 


and  that  at  the  bottom  (z  = — h)  being 


A^-  = kub 
dz 


(3.35) 


where  u is  velocity,  ( is  surface  elevation,  A is  the  vertical  eddy  viscosity,  tw  is  wind 
stress,  p is  water  density,  k is  the  linearized  bottom  friction  coefficient,  and  ub  is  the 
bottom  velocity. 

Integrating  Eq.  (3.32)  twice  and  applying  the  boundary  conditions,  the  solution  is 
given  as  follow  with  a constant  vertical  eddy  viscosity  (Neumann  and  Pierson,  1969) 


(3.36) 


d(.(z2-h2)  h,  tw.  A, 

where  dC^jdx  is  calculated  from  Equations  (3.32)  and  (3.36)  as 

(K  _ 3 tw  ( hk  + 2A ) 
dx  2 pgh  ( hk  + ZA) 

The  following  values  were  used  in  the  model  simulation  for  the  steady  state  wind 
setup 


(3.37) 


tw  — 1000  dyne /cm2 


P = 
k = 
l = 

(H0  + Hi)/ 2 = 

0 < s < 0.1 

At  = 


1 gm/cm3 
1 cm/s 
62  km 
4000  m 

100  s 


(3.38) 


Comparison  between  model  results  and  the  analytic  solution  for  the  setup  as  a 
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HI  NO  SETUP  TEST  HITH  SLOPING  BOTTOM 
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Figure  3.14:  Comparison  between  simulated  maximum  surface  elevation  and  analytic 
solution  for  a sloping  bottom  basin  forced  by  wind  ( a0  is  the  setup  in  the  flat  bottom). 
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Figure  3.15:  Comparison  between  simulated  vertical  velocity  structure  and  analytic 
solution  for  a constant  wind  in  a closed  basin. 
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Figure  3.16:  Grid  system  for  wind  setup  test  in  a bowtie-shaped  basin. 
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Figure  3.17:  Steady  state  setup  in  cm  for  (a)  1 dyne/cm 2 eastward  wind,  and  (b)  1 
dyne/ cm?  northward  wind. 
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function  of  bottom  slope  is  shown  in  Figure  3.14.  The  model  estimates  the  setup  at 
all  slopes  with  less  than  3%  error. 

The  following  values  were  used  in  the  model  simulation  for  the  vertical  structure 
of  horizontal  velocity 


tw 

P 

k 

l 

{H0  + H,)/2 
A 
A t 


1 dyne/ cm2 
1 gm/cm 3 
0.5  cm/s 
62  km 
5 m 

5 cm2 /s 
300  s 


(3.39) 


Comparison  between  model  results  and  analytic  solution  of  the  vertical  structure 
of  horizontal  velocity  is  shown  in  Figure  3.15.  Eleven  vertical  levels  were  used  in  this 
test.  The  error  is  less  than  2 %.  When  the  number  of  vertical  levels  was  increased,  the 
model  approached  the  analytic  solution,  particularly  near  the  bottom.  With  fewer 
vertical  levels,  the  error  was  found  to  increase  especially  near  the  bottom  although 
the  model  results  still  compare  well  with  the  analytic  vertical  structure  of  velocity. 

The  third  test  for  the  steady-state  response  to  wind  forcing  in  a bowtie-shaped 
estuary  was  performed.  Figure  3.16  shows  the  test  grid  system  in  the  bowtie-shaped 
basin.  The  basin  has  a longest  distance  of  60  km  in  the  y-direction,  and  a shortest 
distance  of  25  km  in  the  x-direction.  A constant  wind  stress  of  1 dyne/ cm2  to  the 
north  and  east  was  used  for  the  simulation.  The  bottom  friction  coefficient  was  set 
to  zero.  As  shown  in  Figure  3.17,  the  simulated  setup  is  in  excellent  agreement  with 
the  analytic  solution  in  both  directions. 


98 

3.5  Density-Driven  Currents 


Fresh  water  runoff  and  intrusion  of  salt  water  from  ocean  combine  to  produce 
density-driven  currents.  For  testing  baroclinic  flow  simulation,  analytic  solutions  de- 
veloped by  Hansen  and  Rattray  (1965)  were  used.  A typical  estuarine  flow  shows 
a two-layer  system:  seaward  flow  at  the  surface  and  up-estuary  flow  at  the  bottom. 
Neglecting  nonlinear  terms,  Coriolis,  horizontal  diffusion,  and  y-direction  terms,  one 
can  obtain  the  following  continuity,  momentum,  and  salinity  equations  (Hansen  and 
Rattray,  1965) 


u 


dS 

dx 


du  dw 
dx  dz 

(3.40) 

d(  dS  A d2u 

! — + gaz—  = Av— 
dx  dx  dz * 

(3.41) 

dS  d dSs  d dS 

dz  ~ dx^Kh  dx^  + dz^<v  dz^ 

(3.42) 

with  an  equation  of  state  for  estuarine  water 


p = Po(l  + <*S)  (3.43) 

where  £ is  surface  elevation,  S is  salinity,  Av  is  vertical  eddy  viscosity,  u is  horizontal 
velocity,  w is  vertical  velocity,  Kh  is  horizontal  eddy  diffusivity,  Kv  is  vertical  eddy 
diffusivity,  a is  a constant,  p0  is  density  of  freshwater,  and  g is  gravitational  acceler- 
ation. Eq.  (3.41)  means  that  the  surface  elevation  gradient  plus  the  density  gradient 
balance  the  vertical  diffusion  of  momentum. 

Boundary  conditons  are  no-slip  at  the  bottom,  net  horizontal  transport  equals 
river  flow,  and  zero  normal  salt  flux  at  the  boundaries.  In  the  central  regime  of  the 
estuary,  Hansen  and  Rattray  (1965)  assumed  that  both  the  vertical  salinity  stratifi- 
cation and  velocity  structure  are  independent  of  longitudinal  positions.  The  salinity 


solutions  are 
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where 


_S 

S0 


V ,,  1 . 1.0  1 X 

1 + + ~ 2]  ~ 2^  "3) 

/ <j)dri  + / [ <j>  dr)'  drj) 

Jo  Jo  Jo 


4 = j(2  - 3i)  + 1 f)  - —■(>!  - 3i f + 2r/4) 


(3.44) 


(3.45) 


(3.46) 


where  ij  = z/ D and  £ = Rx/BDI<ho  are  dimensionless  vertical  and  horizontal  co- 
ordinates, respectively,  and  zero  subscripts  indicate  values  at  x — 0.  D is  channel 
depth,  B channel  width,  v is  a constant,  </>  stream  function,  Ra  = gaS0D3 / AvKho 
estuarine  Rayleigh  number,  M = KvKh0B2 / R2  a ratio  of  tidal  mixing  to  the  river 
flow,  and  R river  discharge  rate,  v is  the  positive  root  of 


1680M(1  -v)  = 32*/  + 76^2  + (3.47) 

4o  6 4o 

The  horizontal  exchange  coefficient,  K^,  must  increase  seaward  at  a rate  equal  to 
the  integral  mean  velocity  or  fresh  water  discharge  velocity 


±K  = JL 

dx  h BD 


uf 


(3.48) 


For  this  test,  the  computer  model  code  was  modified  to  incorporate  the  above 
non-uniform  horizontal  eddy  diffusivity,  while  also  using  the  following  parameters 


Kv  = 5 cm2 / s 
Kho  = 105  cm2/s 
Av  = 5 cm2  / s 
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Figure  3.18:  Comparison  between  simulated  results  and  analytic  solution  for  the 
vertical  structure  of  horizontal  velocity  induced  by  density  gradient. 
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a = 6 • 10-4 

So  = 30  ppt  (3.49) 

B = 4 km 
D = 10  m 
R = 2 • 107  cm3/s 
p0  = 1.0  g/ cm3 

From  the  above  information,  v has  a value  of  0.2122.  The  simulation  was  started 
with  an  initial  salinity  distribution  which  has  a horizontal  gradient  of  0.318  ppt /km 
and  the  river  inflow  given  above.  As  shown  in  Figure  3.18,  simulation  velocity  struc- 
ture obtained  with  11  vertical  levels  agree  well  with  the  analytic  solution  with  less 
than  2%. 

3.6  Summary 

Tests  performed  in  this  chapter  indicate  the  basic  soundness  of  the  numerical 
schemes  used  for  many  terms  in  the  three-dimensional  model.  Moreover,  the  effect 
of  a curvilinear  grid  on  accuracy  of  the  model  was  also  tested.  Table  3.1  summarizes 
all  the  test  runs  made  here.  However,  good  agreement  between  model  results  and 
analytic  solutions  of  idealized  linearized  problems  only  ensures  the  numerical  accuracy 
of  the  model  and  does  not  ensure  the  correctness  of  the  parameterization  of  physical 
processes  ( e.g .,  turbulent  mixing,  bottom  friction,  and  nonlinear  acceleration,  etc.). 
To  test  the  soundness  of  the  turbulence  model  and  bottom  friction  formulation,  it  is 
necessary  to  compare  model  results  with  laboratory  and  field  data. 
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Table  3.1:  Characteristics  of  model  analytical  tests. 


Model 

Features 

Forced 
Constant 
Surface  Slope 

Seiche 

Tidal 

Forcing 

Wind 

Forcing 

Density- 

Driven 

Basin 

Geometry 

uniform  and 
non-uniform 
rectangular 

uniform 

rectangular 

uniform 

rectangular, 

annular 

uniform 

rectangular, 

bowtie 

uniform 

rectangular 

Basin 

Bottom 

flat 

flat 

flat, 

sloping 

flat, 

sloping 

flat 

Dimension- 

ality 

2-D 

2-D 

2-D 

2-D, 3-D 

3-D 

Surface 

Slope 

yes 

yes 

yes 

yes 

yes 

Local  Ac- 
celeration 

no 

yes 

yes 

yes 

yes 

Bottom 

Stress 

yes 

no 

ye£ 

yes 

no 

Vertical 

Diffusion 

no 

no 

no 

yes 

yes 

CHAPTER  4 

ANALYSIS  OF  OBSERVED  HYDRODYNAMIC  DATA  FROM  JAMES  RIVER 


Before  the  hydrodynamic  model  is  applied  to  James  River,  it  is  necessary  to  learn 
about  the  circulation  of  James  River  from  observed  hydrodynamic  data.  A compre- 
hensive set  of  hydrodynamic  data  has  been  collected  from  the  James  River  estuary 
by  the  Virginia  Institute  of  Marine  Science  (VIMS).  Byrne  et  al.  (1987)  and  Hep- 
worth  and  Kuo  (1989)  presented  a detailed  description  of  data  collection  techniques 
and  analyses  of  data.  The  major  elements  of  the  data  collection  program  included 
a network  of  tide  guages,  a series  of  current  meters  and  CTD  mooring  deployments 
(Figure  4.1),  and  a series  of  slack  water  measurements  of  vertical-longitudinal  salinity 
profiles.  Time  series  of  measured  surface  elevations  and  currents  were  analyzed  using 
a least-squares  method  to  obtain  the  harmonic  constants.  These  results  will  be  used 
for  model  calibration  and  verification  by  comparing  them  with  model  results.  A series 
of  slack  water  vertical-longitudinal  salinity  profiles  along  the  axis  of  the  river  will  be 
used  to  calibrate  the  vertical  mixing  parameterization  in  the  salinity  equation.  The 
Eulerian  residual  current  distribution  over  the  transect  ‘JB’  shown  in  Figure  4.2  will 
be  used  to  evaluate  the  ability  of  the  model  to  reproduce  the  general  features  of  the 
residual  current  structure  over  this  section. 

In  the  next  chapter,  the  three-dimensional  hydrodynamic  model  will  be  applied 
to  James  River  and  model  results  will  be  compared  with  field  data  described  in  this 
chapter  with  various  sensitivity  tests  to  quantify  the  parameter-varying  sensitivity. 
And  stratification-destratification  phenomena  of  vertical  salinity  structure  will  be 
simulated  according  to  the  spring-neap  tidal  cycle.  Next  salinity  intrusion  problem 
and  frontal  system  will  be  analyzed.  And  model  calibration  and  validation  process 
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Table  4.1:  Tidal  amplitude  harmonic  constants  obained 
from  June  19  to  July  17,  1985.  (period  in  hour,  ampli- 
tude in  cm,  phase  lag  in  degree) 


Station  No. 

1 

2 

3 

4 

NOS  No. 

Const. 

Period 

1 

M2 

12.42 

Amp. 

33.4 

36.4 

36.7 

37.0 

Phase  lag 

327.6 

313.9 

286.9 

278.9 

2 

S2 

12.00 

Amp. 

4.4 

4.5 

5.0 

5.0 

Phase  lag 

338.3 

318.1 

278.2 

267.8 

3 

N2 

12.66 

Amp. 

5.9 

6.3 

6.7 

6.6 

Phase  lag 

193.0 

185.2 

152.0 

143.9 

4 

K1 

23.93 

Amp. 

4.4 

4.3 

4.4 

4.5 

Phase  lag 

327.0 

310.9 

296.3 

291.4 

5 

M4 

6.21 

Amp. 

1.0 

0.8 

0.5 

0.5 

Phase  lag 

183.2 

162.4 

83.2 

310.0 

6 

01 

25.82 

Amp. 

4.6 

4.5 

4.9 

4.9 

Phase  lag 

22.5 

5.9 

347.1 

341.7 

7 

M6 

4.14 

Amp. 

1.0 

0.2 

0.4 

0.5 

Phase  lag 

18.1 

359.5 

0.3 

304.2 

wil  be  discussed.  Finally  we  will  discuss  the  a-  and  2-baroclinic  effects. 

4.1  Water  Level  Data 

According  to  the  field  data  analysis,  the  sea  surface  elevation  records  showed  that 
tidal  forcing  accounts  for  about  85  to  95  % of  the  surface  fluctuations  in  the  James 
River.  Clearly,  tide  is  the  dominant  hydrodynamic  factor.  The  semi-diurnal  M2 
constituent  was  the  strongest  component  with  the  form  ratio,  F = (K 1 +Ol)/(M2  + 
S 2)  = 0.2,  which  is  classified  as  a semi-diurnal  dominant  tidal  type.  The  other 
constituents  contribute  only  variations  such  as  the  spring-neap  cycle.  Shallow  water 
components  of  M4  and  M6,  caused  from  M2,  were  negligible  with  less  than  1 cm  in 
most  areas.  The  tidal  harmonic  analysis  for  a 29-day  period  starting  from  June  19, 
1985  is  shown  in  Table  4.1.  The  M2  constituent  at  the  open  boundary  of  Station  4 
has  an  amplitude  of  37.0  cm,  and  is  modulated  by  the  S2  (5.0  cm)  and  N2  (6.6  cm) 
constituents  at  periods  of  14.8  and  27.6  days,  respectively.  The  interaction  of  the 
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Figure  4.1:  Lower  James  River  and  Hampton  Roads  (After  Hepworth  and  Kuo 
(1989)). 


Depth  (in  meters) 
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Figure  4.2:  Observed  Eulerian  residual  velocity  over  transect  ‘JB’  (After  Hamrick  et 
al.  ,1989). 
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three  constituents  causes  the  semi-diurnal  tidal  amplitude  to  vary  between  25.4  cm 
at  neap  condition  to  48.6  cm  at  spring  condition. 

The  amplitudes  of  the  two  significant  diurnal  constituents  01  (4.9cm)  and  K1 
(4.5  cm)  are  about  25  % of  the  M2  constituent.  Their  only  significant  contribution 
to  the  tidal  character  of  the  area  is  to  create  a diurnal  inequality,  which  is  shown  in 
Figure  4.3. 

4.2  Tidal  Currents 

Tidal  currents  in  the  James  River  are  again  dominated  by  the  M2  constituent, 
which  have  a surface  amplitude  of  50.5  cm/s  in  u-  direction  and  38.8  cm/s  in  in- 
direction at  Station  A of  ‘JB’  in  Figure  4.1.  The  near  bottom  currents  are  half  as 
high  as  the  surface  currents.  The  S2  and  N2  constituents  have  approximately  one  sixth 
the  magnitude  of  the  M2  and  interact  with  it  to  create  the  spring-neap  variation  in 
the  semidiurnal  band.  The  magnitude  of  the  diurnal  constituents  (Kl,  01)  are  quite 
significant  on  the  order  of  7 - 12  cm/s,  while  the  M4  and  M6  are  quite  small,  on  the 
order  of  1 - 2 cm/s. 

4.3  Circulation  Pattern 

Circulation  patterns  in  the  James  River  estuary  were  reported  based  on  hydraulic 
model  experiments  (Byrne  et  al.,  1987).  They  showed  that  the  lower  James  River 
was  a “leaky  closed  system”  wherein  surface  circulation  is  expanded  as  an  elongated 
counterclockwise  gyre  with  slight  losses  of  material  upstream  past  Mulberry  Point 
and  substantial  losses  downstream  past  Sewells  Point.  Bottom  circulation  has  a net 
upstream  progression  with  some  indication  of  counterclockwise  tendency  in  the  upper 
reaches  of  Burwell  Bay.  Nontidal  circulation  pattern  was  substantiated  by  results  of 
a current  measurement  experiment  conducted  near  the  James  River  Bridge  in  the 
summer  of  1985  (Byrne  et  al.,  1987).  The  longitudinal  component  of  residual  currents 
shows  a circulation  at  this  location  as  a typical  two  layer  system  (Figure  4.2).  The 
level  of  no  net  motion  in  Figure  4.2  is  higher  on  the  left  when  looking  downstream. 
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Circulation  in  Hampton  Roads  based  on  hydraulic  model  experiments  indicates  the 
cyclonic  nature  of  circulation. 

4.4  Slack-Water  Salinity  Data 

A series  of  slack  water  vertical-longitudinal  salinity  profiles  along  the  axis  of  the 
river  is  shown  in  Figure  4.4  - 4.7.  These  will  be  used  to  calibrate  the  vertical  mixing 
parameterization  in  the  salinity  equation. 

4.5  Front  Dynamics 

Fronts  are  regions  where  flow  convergence  and  strong  exchange  between  different 
water  masses  occurs,  and  are  hence  important  in  estuarine  dynamics.  Estuarine  fronts 
have  significant  biological  consequences,  since  they  are  areas  of  high  productivity  for 
all  elements  of  the  food  chain  from  phytoplankton  to  fish.  The  location  of  estuarine 
frontal  zones  should  be  considered  in  the  design  of  waste  water  discharge  outfalls 
to  avoid  excessive  concentration  and  contamination.  Frontal  circulation  affects  the 
dispersion  of  oil  spills  and  ocean  dumped  wastes.  Thus,  the  study  of  estuarine  fronts 
must  be  incorporated  into  environmental  monitoring  programs  to  determine  how  pol- 
lutants are  transported  and  concentrated  into  biologically  active  areas. 

Frontogenesis  in  estuaries  is  often  a transient  and  local  phenomenon  which  is 
related  to  particular  geometric  features  and  occurring  during  a particular  portion  of 
a tidal  cycle.  In  coastal  areas  where  tidal  currents  are  strong  enough,  outflow  plumes 
may  be  swept  back  into  the  estuary  by  the  flood  current.  Simpson  and  Nunes  (1981) 
called  it  a tidal  intrusion  front.  According  to  them,  the  frontal  boundary,  marking 
the  surface  convergence  where  the  incoming  water  flows  under  the  less  dense  surface 
layer,  propagates  upriver  in  a characteristic  V-shaped  configuration  with  strong  point 
convergence  at  the  apex  of  the  V. 

Observations  revealed  that  the  front  in  the  lower  James/Hampton  Roads  is  V- 
shaped  and  persistent  in  terms  of  its  occurrence  (Byrne,  1987).  The  front  occurs  in  a 
small  area  off  Newport  News  Point  at  flood  tide  during  every  tidal  cycle,  although  its 
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ELEVATION  AT  OPEN  BOUNDARY  (JUNE  19  - JULY  17.1985) 


TIME  (HOURS) 


Figure  4.3:  Surface  elevation  at  the  open  boundary  between  June  19  and  July  17, 
1985. 
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Figure  4.4:  Observed  longitudinal- vertical  salinity  at  slackwater  before  flood  on  June 
19,  1985. 
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Figure  4.5:  Observed  longitudinal-vertical  salinity  at  slackwater  before  flood  on  July 
3,  1985. 
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Figure  4.6:  Observed  longitudinal-vertical  salinity  at  slackwater  before  flood  on  July 
9,  1985. 
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Figure  4.7:  Observed  longitudinal- vertical  salinity  at  slackwater  before  flood  on  July 
17,  1985. 
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Figure  4.8:  Configuration  of  fronts  at  1118  hours.  From  aerial  photographs  of  August 
29,  1985.  Generalized  flow  direction  shown  by  arrows.  Dotted  lines  are  depth  contours 
in  ft.  (After  Byrne,  1987) 
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Figure  4.9:  Configuration  of  fronts  at  1400  hours.  From  aerial  photographs  of  August 
29,  1985.  Generalized  flow  direction  shown  by  arrows.  Dotted  lines  are  depth  contours 
in  ft.  (After  Byrne,  1987) 
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evolution  with  time  is  somewhat  variable  with  respect  to  tide  condition.  During  its 
life  span  the  front  is  revealed  as  a narrow  zone  of  surface  flow  convergence  and  strong 
vertical  transport.  The  vertical  transport  across  the  frontal  system  was  examined 
by  the  dye  release  experiment  (Byrne,  1987).  The  dye  was  released  at  the  surface 
near  the  downriver  side.  When  the  upriver  water  column  was  stratified,  the  peak  dye 
concentration  tended  to  locate  at  depth  where  the  vertical  gradient  was  greatest.  This 
shows  that  the  material  in  the  upper  portion  of  the  water  column  are  transported 
into  lower  horizons  where  the  net  estuarine  transport  is  upriver.  This  also  tells  that 
there  is  a strong  vertical  flow. 

According  to  Byrne  (1987),  field  observations  suggest  that  the  frontal  system  and 
its  general  evolution  may  be  described  in  terms  of  three  segments  as  shown  in  Fig- 
ure 4.8.  The  initial  phase  in  frontal  development  is  associated  with  flow  exiting  Hamp- 
ton Flats  near  Newport  News  Point,  designated  as  branch  la  of  the  frontal  system. 
During  the  early  phase  of  flood  tide,  there  is  also  a distinct  boundary  parabathytic 
with  the  southern  flank  of  Hampton  Flats.  This  boundary  designated  as  branch  lb, 
apparently  represents  flow  contributions  from  flow  coursing  between  Hampton  Flats 
and  the  navigation  channel.  These  two  segments  usually  join  together  at  the  western 
edge  of  Newport  News  Bar  to  form  a V shape  front.  The  front  position  generally  rests 
near  the  9 m depth  contour.  The  third  segment,  branch  II,  of  the  frontal  system  is 
that  associated  with  the  advancing  flood  current  along  the  Newport  News  Channel. 
So  as  the  flood  tide  currents  progress  up  the  channel  section,  the  more  saline  waters 
entering  Hampton  Roads  encounter  the  estuarine  salinity  gradient.  Another  bound- 
ary segment  is  formed  which  progresses  upriver  with  the  flood  tide,  and  frequently 
joins  an  apex  formed  by  the  three  boundaries  on  the  slope  leading  to  the  deep  water 
of  the  lower  James  (Figure  4.9). 


CHAPTER  5 

APPLICATION  OF  THREE-DIMENSIONAL  MODEL  TO  JAMES  RIVER 


5.1  Procedure  of  Model  Simulation 

In  this  chapter,  model  application  to  James  River  is  described.  First  of  all,  the 
model  simulation  procedure  is  described.  Simulation  of  tidal  circulation  is  then  pre- 
sented, followed  by  simulation  of  the  stratification- destratification  process  cycle  as- 
sociated with  the  spring-to-neap  cycle.  Following  that,  simulation  of  salt  intrusion  is 
presented  before  a detailed  presentation  on  the  simulation  of  formation  and  evaulation 
of  a tidal  intrusion  front. 

First  of  all,  using  the  grid  generation  program,  two  sets  of  non-orthogonal  grid 
systems  in  the  horizontal  plane  were  generated  for  the  James  River:  a coarse  grid 
system  with  76x38  grid  points  (Figure  5.1)  and  a fine  grid  system  with  121x45  (Fig- 
ure 5.2).  Figure  5.3  and  Figure  5.4  show  the  same  grid  systems  when  plotted  in  (£,r/) 
coordinates,  respectively.  A space-staggered  grid  system  as  discussed  before  was  used. 
Next,  the  depth  at  each  grid  point  was  interpolated  from  bathymetric  data  on  the 
NOAA/NOS  navigation  chart.  In  the  vertical  direction,  a total  of  11  vertical  levels 
is  used  with  a uniformly  spaced  a grid. 

The  simulations  performed  here  use  a two-time-level  scheme  for  time  stepping, 
and  generally  central  difference  schemes  for  spatial  gradients.  The  external  mode  is 
solved  implicitly,  and  the  internal  mode  is  solved  explicitly  except  that  the  vertical 
diffusion  term  is  treated  implicitly.  The  vertical  eddy  coefficients  for  momentum  and 
salinity  are  determined  by  the  simplified  second-order  closure  schemes  of  Sheng  (1983) 
and  Sheng  and  Villaret  (1989). 

The  model  was  applied  to  simulate  the  circulation  and  salinity  transport  in  the 
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James  River  during  a 29-day  period  starting  June  19,  1985.  Based  on  a series  of 
four  synoptic  scale  slack-water  (slack-before-flood  (SBF))  salinity  surveys  made  in 
the  James  River  between  June  19  and  July  17,  1985,  it  is  apparent  that  the  estuary 
goes  through  a periodic  cycle  of  stratification  and  destratification  according  to  the 
neap-spring  tidal  variation.  Salinity  was  measured  at  slack-before-flood,  from  the 
mouth  of  the  river  to  the  head  of  the  salt  intrusion  with  5 to  10  km  distance  between 
two  neighboring  stations.  Longitudinal  and  vertical  salinity  distributions  for  each 
survey  are  shown  in  chapter  4. 

To  simulate  the  process  of  stratification  and  destratification,  it  is  necessary  to 
define  river  discharge,  tidal  forcing  at  the  open  boundary,  and  initial  conditions  of 
salinity  and  flow  field.  The  river  discharge  at  Richmond  averaged  81  m3/sec  during 
the  sampling  period,  and  this  value  was  used  for  the  fresh  water  flow  in  the  model 
simulation.  For  tidal  forcing  at  the  open  boundary  ( i.e .,  entrance  to  Chesapeake 
Bay),  7 tidal  harmonic  constituents  (NOS  number  : 1 to  7)  were  used  and  computed 
by  the  method  of  least  squares  (Boon  and  Kiley,  1978).  A time  series  of  surface  level 
at  the  open  boundary  is  shown  in  Figure  4.3.  For  salinity  boundary  conditions,  a 
constant,  vertically- uniform  salinity  of  26  ppt  was  used  during  the  simulation. 

For  the  initial  conditions  of  salinity  and  flow  field,  the  model  is  first  spun  up  by  a 
M2  tide  with  an  amplitude  of  38  cm  at  the  open  boundary  which  corresponds  to  the 
tidal  condition  on  June  19,  1985.  The  model  is  spun  up  from  the  condition  of  no  flow 
and  homogeneous  temperature  for  5 days  without  considering  salinity.  The  model 
simulation  is  then  continued  for  5 more  days  with  an  initial  salinity  field  interpolated 
from  the  slackwater  data  on  June  19,  1985.  During  this  5-day  simulation,  the  salinity 
field  is  held  fixed  while  the  velocity  field  is  allowed  to  adjust  to  the  baroclinic  field. 
After  the  5-day  baroclinic  adjustment,  the  model  is  allowed  to  start  the  dynamic 
simulation  from  June  19  to  July  17,  1985,  wherein  both  the  velocity  and  salinity  field 
can  change  with  time. 
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Figure  5.1:  Curvilinear  grid  for  the  James  River  coarse  grid  system  (physical  domain). 


120 


Figure  5.2:  Curvilinear  grid  for  the  James  River  fine  grid  system  (physical  domain). 
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Figure  5.3:  Curvilinear  grid  for  the  James  River  coarse  grid  system  (computation 
domain). 
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Figure  5.4:  Curvilinear  grid  for  the  James  River  fine  grid  system  (computation  do- 
main). 
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5.2  Simulation  of  Tidal  Circulation 


Comparison  between  simulated  (coarse  grid)  and  tidal  elevation  and  tidal  currents 
was  made  by  means  of  harmonic  analysis  of  the  computed  and  observed  data.  The 
computed  and  observed  water  levels  at  the  3 stations  shown  in  Figure  5.1  were  ana- 
lyzed using  the  harmonic  analysis  program  and  the  results  were  compared  with  data. 
Table  5.1  lists  the  amplitudes  and  phase  angles  for  5 tidal  constituents  at  Stations  1, 
2,  and  3.  It  is  apparent  that  the  computed  amplitude  and  phase  angle  agree  well  with 
observations  with  a maximum  error  of  9 % in  surface  elevation  amplitude  at  Station 
3 and  with  a maximum  error  of  6 minutes  of  phase  lag  at  Station  2 (M2  case). 

The  harmonic  constants  of  computed  tidal  currents  were  also  compared  with  the 
observed  ones.  Table  5.2  lists  the  amplitude  and  phase  angle  of  computed  tidal  cur- 
rents at  four  different  levels  at  Station  4.  The  model  generally  predicts  the  amplitude 
of  near-surface  currents  quite  well  within  2 % error  (M2  case),  but  appears  to  over- 
predict the  amplitude  of  near-bottom  currents.  The  error  may  be  due  to  the  use  of 
a relatively  poorly  resolved  bottom  topography  associated  with  the  coarse  grid  and 
due  to  an  inaccurate  density-driven  circulation. 
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Figure  5.5:  Eulerian  residual  current  at  surface  in  the  James  River  between  June  19 
and  July  17,  1985. 
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Figure  5.6:  Eulerian  residual  current  at  mid-depth  in  the  James  River  between  June 
19  and  July  17,  1985. 
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Figure  5.7:  Eulerian  residual  current  near  bottom  in  the  James  River  between  June 
19  and  July  17,  1985. 
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Figure  5.8:  Vertical  velocity  profile  at  spring  near  the  James  River  Bridge.  The 
numbers  on  top  of  plot  indicate  the  time  evolution  in  hours. 
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Figure  5.9:  Vertical  velocity  profile  at  neap  near  the  James  River  Bridge.  The  numbers 
on  top  of  plot  indicate  the  time  evolution  in  hours. 
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TIDAL  CURRENT  NEAR  JAMES  RIVER  BRIDGE 


U-VELOCITY  (CM/SEC) 


Figure  5.10:  Computed  near-surface  tidal  currents  near  the  James  River  Bridge  be- 
tween June  19  and  July  17,  1985. 
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Figure  5.11:  Contours  of  29-day  residual  longitudinal  velocity  at  a transect  near  the 
James  River  Bridge  (along  i=65). 
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Table  5.1:  Comparison  of  the  model  predicted  harmonic  constants  of  water  level  with 
the  measured  data. 


Constituent 

M2 

S2 

N2 

K1 

01 

St.l 

0=57, 

j=19) 

Amplitude 

(cm) 

Data 

33.4 

4.4 

5.9 

4.4 

4.6 

Model 

36.4 

5.6 

6.2 

2.9 

4.5 

Phase  lag 
(degree) 

Data 

327.6 

338.3 

193.0 

327.0 

22.5 

Model 

327.0 

315.5 

190.1 

316.5 

7.8 

St.2 

0=64, 

j=16) 

Amplitude 

(cm) 

Data 

36.4 

4.5 

6.3 

4.3 

4.5 

Model 

37.5 

6.0 

6.5 

2.9 

4.6 

Phase  lag 
(degree) 

Data 

310.4 

299.4 

174.6 

308.4 

359.0 

Model 

313.4 

303.9 

179.0 

311.4 

2.9 

St.3 

(i=74, 

j=17) 

Amplitude 

(cm) 

Data 

36.7 

5.0 

6.7 

4.4 

4.9 

Model 

36.5 

5.9 

6.4 

2.8 

4.5 

Phase  lag 
(degree) 

Data 

286.9 

278.2 

152.0 

296.3 

347.1 

Model 

286.1 

265.5 

150.3 

294.2 

347.4 

The  simulated  Eulerian  residual  currents  over  the  entire  James  River  at  surface- 
level,  mid-depth,  and  bottom-level  during  the  29-day  simulation  are  shown  in  Fig- 
ures 5.5,  5.6,  and  5.7,  respectively.  The  results  exhibit  a typical  estuarine  circulation 
pattern,  showing  a down-estuary  flow  in  the  upper  portion  of  the  water  column  and 
an  up-estuary  flow  in  the  lower  portion  of  the  water  column.  Generally  the  velocity 
differences  through  the  water  column  associated  with  an  estuarine  circulation  will  be 
an  order  of  magnitude  greater  than  the  velocity  differences  associated  with  a riverine 
shear  flow  in  which  there  is  no  longitudinal  density  gradient  force  (Officer,  1978). 
The  shear  in  the  vertical  velocity  profile  is  due  to  a combination  of  bottom  friction 
and  internal  friction  effects.  In  the  middle  of  lower  reaches  of  most  estuaries,  how- 
ever, the  vertical  velocity  profile  is  due  to  a combination  of  the  longitudinal  surface 
slope  (barotropic)  term  and  the  longitudinal  density  gradient  (baroclinic)  term  which 
acts  in  the  opposite  direction  and  increases  as  a function  of  depth.  Figure  5.8  and 
Figure  5.9  show  the  computed  vertical  structure  of  the  horizontal  velocities  near  the 
James  River  Bridge  (Station  4)  at  spring  and  neap  tides,  respectively.  Two-layer 
estuarine  flows  are  clearly  seen  in  these  plots  with  more  positive  u-velocities  in  the 
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Table  5.2:  Comparison  of  the  model  predicted  harmonic  constants  of  velocity  with 
the  measured  data  (St.  4,  Depth=8.5  m) 


depth 

Const. 

M2 

S2 

N2 

K1 

01 

1.2  m 

Amplitude 

(cm/sec) 

u 

Data 

50.5 

7.2 

9.8 

10.1 

6.0 

Model 

49.5 

7.2 

8.0 

4.5 

5.7 

V 

Data 

38.8 

4.2 

5.7 

8.4 

4.3 

Model 

40.9 

5.4 

5.8 

5.0 

5.9 

Phase  lag 
(degree) 

u 

Data 

128.4 

147.3 

196.5 

106.7 

167.7 

Model 

120.2 

117.8 

188.0 

77.4 

144.8 

V 

Data 

308.0 

327.0 

18.1 

282.2 

343.6 

Model 

323.2 

314.9 

22.0 

293.6 

353.2 

3.1  m 

Amplitude 

(cm/sec) 

u 

Data 

48.8 

7.0 

7.3 

8.5 

5.9 

Model 

48.1 

6.7 

7.6 

4.5 

1.9 

V 

Data 

28.6 

4.7 

4.1 

4.4 

2.4 

Model 

36.7 

4.8 

5.5 

4.0 

3.7 

Phase  lag 
(degree) 

u 

Data 

123.5 

143.3 

168.8 

106.7 

168.2 

Model 

125.8 

120.4 

188.7 

75.9 

140.7 

V 

Data 

294.9 

319.2 

355.2 

263.4 

324.8 

Model 

316.3 

303.7 

10.7 

256.7 

321.7 

6.3  m 

Amplitude 

(cm/sec) 

u 

Data 

38.6 

4.4 

6.1 

5.8 

4.5 

Model 

34.2 

4.5 

4.9 

2.5 

2.9 

V 

Data 

10.2 

1.3 

2.2 

1.7 

1.7 

Model 

24.9 

2.9 

2.9 

1.0 

1.5 

Phase  lag 
(degree) 

u 

Data 

112.4 

137.5 

182.2 

94.4 

149.3 

Model 

113.3 

94.6 

162.6 

62.7 

131.6 

V 

Data 

283.2 

329.8 

340.5 

275.9 

353.3 

Model 

261.8 

295.5 

329.7 

223.6 

284.6 

8.5  m 

Amplitude 

(cm/sec) 

u 

Data 

30.3 

3.2 

5.2 

4.1 

3.1 

Model 

24.3 

3.6 

3.8 

2.1 

2.1 

V 

Data 

6.5 

1.0 

1.8 

1.2 

1.1 

Model 

16.2 

2.3 

2.3 

0.7 

0.9 

Phase  lag 
(degree) 

u 

Data 

105.8 

136.5 

179.6 

95.5 

151.1 

Model 

111.91 

95.8 

164.3 

63.8 

135.3 

V 

Data 

285.1 

326.6 

357.3 

289.9 

310.4 

Model 

256.9 

258.6 

327.7 

230.9 

293.5 
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upper  layer  and  more  negative  u-velocities  in  the  lower  layer.  During  spring  tide, 
velocity  profiles  show  a larger  range  of  fluctuation  due  to  the  high  tidal  range. 

Instantaneous  tidal  currents  at  station  4 near  the  James  River  Bridge  during  the 
29-day  period  are  shown  in  Figure  5.10,  which  agree  quite  well  with  the  data  shown 
in  Hamrick  et  al.  (1990). 

An  isotach  contour  plot  of  simulated  residual  current  across  the  James  River 
Bridge  transect  are  shown  in  Figure  5.11.  The  isotach  lines  compared  quite  well  with 
the  measured  data,  considering  the  fact  that  there  are  only  four  grid  points  across 
the  James  River  in  this  section.  The  residual  isotach  lines  are  very  sensitive  to  the 
bathymetry  and  horizontal  boundary  conditions  of  flow  and  salinity. 

According  to  Pritchard  (1956)  and  Hamrick  (1989),  the  net  non-tidal  velocity 
showed  a down-estuary  movement  in  the  upper  half  of  the  water  column  and  an 
up-estuary  movement  in  the  lower  half  of  the  water  column.  Figure  4.2  shows  the 
two-layer  circulation  pattern  near  the  James  River  Bridge,  where  the  mean  vertical 
salinity  variation  was  about  4 ppt  in  9 m depth.  The  net  circulation  in  James  River  is 
primarily  due  to  the  freshwater  discharge  and  the  gravitational  convection  caused  by 
the  density  difference  between  freshwater  and  seawater.  The  gravitational  convection 
is  a result  of  dynamic  interactions  between  the  salinity  and  current  distributions,  and 
tends  to  increase  stratification. 

5.3  Sensitivity  Tests 

The  purpose  of  sensitivity  tests  is  to  determine  the  sensitivity  of  the  results  of  the 
present  curvilinear  hydrodynamic  modeling  system  to  changes  in  model  parameters, 
and  thereby  to  gain  insight  into  the  relevance  of  the  model  parameters  to  the  physical 
processes  governing  salinity  transport. 

A total  of  12  test  cases  were  run  with  various  combinations  of  model  parameters 
shown  in  Table  5.3.  The  sensitivity  of  model  results  to  changes  in  the  various  model 
parameters  was  investigated  by  comparing  results  of  the  test  cases  (Table  5.3)  with 


134 

measured  data  or  the  baseline  run  (case2)  which  was  presented  in  the  last  section. 
Table  5.5,  Table  5.6,  and  Table  5.7  list  the  tidal  elevation  constituents  at  stations  1, 
2,  and  3,  respectively.  Table  5.8  lists  the  tidal  constituents  of  near-surface  u-velocity 
(east- west  direction)  at  station  4.  Table  5.9  lists  the  tidal  constituents  of  near-surface 
u-velocity  at  station  4. 

5.3.1  Effect  of  Bottom  Friction 

Bottom  friction  is  frequently  the  most  important  non-linear  term  in  shallow  water 
and  generates  overtides  on  its  own  account  (Doodson  and  Proudman,  1923).  Bottom 
friction  influences  both  the  tidal  range  and  the  phase  lag.  It  is  possible  to  adjust 
the  bottom  friction  coefficient  to  achieve  better  agreement  between  predicted  and 
measured  surface  elevations.  The  zone  where  the  flow  is  appreciably  retarded  by 
friction  against  the  surface  and/or  bottom  is  known  as  the  boundary  layer.  The  most 
obvious  constraint  is  that  the  viscosity  of  the  fluid  enforces  the  no-slip  condition. 
Bottom  stress  can  be  expressed  as  a linear  or  quadratic  function  of  the  velocity. 
Normally  in  tidal  dynamics  a quadratic  law  is  assumed.  Bottom  stress  is  obtained  by 
matching  the  measured  velocity  profile  to  the  logarithmic  law  of  the  wall  in  shallow 
water.  The  logarithmic  velocity  profile  in  the  inertial  sublayer  is  one  of  the  major 
results  in  turbulent  theory.  From  the  law  of  the  wall,  the  drag  coefficient  Cj  is 
determined  by 

Cd  = n2[\n(H  + zb)/z0]~2  (5.1) 

where  zb  is  the  grid  point  nearest  the  bottom,  z0  is  the  roughness  height,  and  k is  the 
von  Karman  constant  set  to  0.40  here. 

Bursinger  et  al.  (1971)  also  defined  a vertically-varying  Monin- Oboukhov  length 
as 

, - u r ,,o, 

ngw'p'/po  K,gl<v(dp/dz)/p0 

where  is  the  frictional  velocity  defined  as  (rj/p)1/2,  where  rb  is  the  bottom  shear 
stress, p is  the  density,  and  Kv  is  the  vertical  eddy  diffusivity.  They  calculated  a 
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Table  5.3:  Test  cases  performed. 


Grid 

system 

Bottom 

friction 

Bottom  rough- 
ness height 

Monin- 

Obukhov 

Law  of 
the  wall 

Constant 
bottom  friction 

easel 

coarse 

explicit 

0.4 

yes 

yes 

no 

case2 

coarse 

implicit 

0.4 

yes 

yes 

no 

case3 

coarse 

implicit 

0.4 

no 

yes 

no 

case4 

coarse 

implicit 

0.1 

yes 

yes 

no 

case5 

coarse 

implicit 

0.04 

yes 

yes 

no 

case6 

coarse 

implicit 

no 

no 

no 

0.002 

case7 

coarse 

implicit 

no 

no 

no 

0.003 

case8 

coarse 

implicit 

no 

no 

no 

0.004 

case9 

fine 

implicit 

0.4 

yes 

yes 

no 

case 10 

fine 

implicit 

0.04 

yes 

yes 

no 

easel 1 

fine 

implicit 

no 

no 

no 

0.003 

casel2 

fine 

implicit 

0.4 

no 

yes 

no 

stability  function  4>(z/ L ) based  on  data  in  a stably  stratified  atmospheric  surface 
layer: 

l+4.7f  (5.3) 

L u*  oz  L 

Monin- Oboukhov  have  successfully  used  z/L  as  the  basic  independent  variable 
for  the  description  of  the  surface  layer,  both  in  stable  and  unstable  conditions.  In 
this  stably  stratified  case,  Cd  is  given  by 

n _ H + zb  , 4.7zn_2 
W — k (In 1 — j 

Zq  L 


(5.4) 
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Table  5.4:  Comparison  between  observed  and  simulated 
tidal  elevation  constituents  at  Station  1 (amplitude  in 
cm,  phase  lag  in  degree). 


constituent 

M2 

S2 

N2 

K1 

01 

data 

amplitude 
phase  lag 

33.4 

4.4 

5.9 

4.4 

4.6 

327.6 

338.3 

193.0 

327.0 

22.5 

easel 

amplitude 
phase  lag 

39.1 

6.1 

6.6 

3.2 

4.9 

328.1 

313.9 

188.9 

310.8 

2.9 

case2 

amplitude 
phase  lag 

33.6 

4.9 

5.4 

2.5 

4.1 

331.1 

323.3 

197.0 

322.5 

13.9 

case3 

amplitude 
phase  lag 

38.3 

5.7 

6.1 

2.9 

4.6 

329.4 

317.4 

191.0 

314.6 

6.7 

case4 

amplitude 
phase  lag 

35.6 

5.4 

5.9 

2.7 

4.4 

328.1 

318.2 

192.1 

318.3 

10.0 

case5 

amplitude 
phase  lag 

36.4 

5.6 

6.2 

2.9 

4.5 

327.0 

315.5 

190.1 

316.5 

7.8 

case6 

amplitude 
phase  lag 

39.1 

6.4 

7.0 

3.4 

5.1 

329.3 

313.1 

188.7 

308.3 

359.5 

case7 

amplitude 
phase  lag 

39.1 

6.3 

6.8 

3.2 

5.0 

328.3 

313.7 

188.6 

309.9 

1.7 

case8 

amplitude 
phase  lag 

38.9 

6.1 

6.5 

3.1 

4.9 

328.1 

313.9 

188.4 

311.7 

3.7 

case9 

amplitude 
phase  lag 

29.3 

3.9 

4.5 

2.1 

3.4 

342.4 

341.3 

212.9 

336.0 

28.4 

case 10 

amplitude 
phase  lag 

32.2 

4.7 

5.3 

2.5 

4.0 

332.7 

328.8 

200.2 

332.5 

23.6 

easel  1 

amplitude 
phase  lag 

32.5 

4.9 

5.5 

2.6 

4.2 

329.9 

324.6 

196.5 

330.5 

21.6 

137 


Table  5.5:  Comparison  between  observed  and  simulated 
tidal  elevation  constituents  at  Station  2 (amplitude  in 
cm,  phase  lag  in  degree). 


constituent 

M2 

S2 

N2 

K1 

01 

data 

amplitude 
phase  lag 

36.4 

4.5 

6.3 

4.3 

4.5 

313.9 

318.1 

185.2 

310.9 

5.9 

easel 

amplitude 
phase  lag 

38.2 

6.3 

6.8 

3.1 

4.8 

310.5 

298.2 

172.8 

305.7 

356.2 

case2 

amplitude 
phase  lag 

35.6 

5.4 

5.8 

2.6 

4.2 

314.3 

305.4 

180.1 

312.3 

3.2 

case3 

amplitude 
phase  lag 

37.7 

6.0 

6.5 

2.9 

4.6 

312.6 

301.4 

175.0 

308.8 

359.9 

case4 

amplitude 
phase  lag 

37.0 

5.8 

6.2 

2.8 

4.5 

311.6 

301.4 

176.5 

309.8 

0.5 

case5 

amplitude 
phase  lag 

37.5 

6.0 

6.5 

2.9 

4.6 

310.4 

299.4 

174.6 

308.4 

359.0 

case6 

amplitude 
phase  lag 

37.8 

6.4 

6.9 

3.2 

4.9 

309.9 

296.7 

171.8 

303.6 

354.2 

case7 

amplitude 
phase  lag 

38.0 

6.4 

6.8 

3.1 

4.8 

310.1 

297.6 

172.4 

304.9 

355.8 

case8 

amplitude 
phase  lag 

37.9 

6.3 

6.7 

3.0 

4.8 

310.6 

298.2 

172.8 

306.0 

356.8 

case9 

amplitude 
phase  lag 

33.9 

4.8 

5.5 

2.4 

3.9 

321.3 

315.1 

190.2 

317.7 

9.0 

case 10 

amplitude 
phase  lag 

36.6 

5.6 

6.2 

2.6 

4.2 

313.5 

307.0 

180.9 

315.7 

6.1 

easel  1 

amplitude 
phase  lag 

36.7 

5.7 

6.2 

2.7 

4.3 

311.1 

304.3 

177.8 

315.0 

5.2 
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Table  5.6:  Comparison  between  observed  and  simulated  tidal  elevation  constituents 
at  Station  3 (amplitude  in  cm,  phase  lag  in  degree). 


constituent 

M2 

S2 

N2 

K1 

01 

data 

amplitude 
phase  lag 

36.7 

5.0 

6.7 

4.4 

4.9 

286.9 

278.2 

152.0 

296.3 

347.1 

easel 

amplitude 
phase  lag 

36.7 

5.9 

6.5 

2.9 

4.5 

285.5 

264.7 

149.4 

293.2 

346.8 

case2 

amplitude 
phase  lag 

35.6 

5.6 

6.2 

2.7 

4.4 

287.7 

267.0 

151.7 

295.2 

348.5 

case3 

amplitude 
phase  lag 

36.3 

5.8 

6.4 

2.8 

4.5 

286.6 

266.2 

150.7 

295.0 

347.6 

case4 

amplitude 
phase  lag 

36.3 

5.8 

6.3 

2.8 

4.4 

286.6 

266.1 

150.9 

294.6 

347.7 

case5 

amplitude 
phase  lag 

36.5 

5.9 

6.4 

2.8 

4.5 

286.1 

265.5 

150.3 

294.2 

347.4 

case6 

amplitude 
phase  lag 

36.9 

6.0 

6.5 

2.9 

4.5 

284.9 

264.1 

148.4 

293.3 

346.0 

case7 

amplitude 
phase  lag 

36.8 

6.0 

6.5 

2.9 

4.5 

285.2 

264.5 

149.1 

293.5 

346.5 

case8 

amplitude 
phase  lag 

36.6 

5.9 

6.5 

2.8 

4.5 

285.6 

264.9 

149.4 

293.7 

346.8 

case9 

amplitude 
phase  lag 

38.2 

6.2 

6.7 

2.8 

4.5 

289.7 

276.9 

153.7 

296.2 

347.7 

case 10 

amplitude 
phase  lag 

38.9 

6.3 

6.8 

2.9 

4.5 

287.8 

267.7 

151.7 

295.3 

347.1 

easel  1 

amplitude 
phase  lag 

38.8 

6.2 

6.8 

2.9 

4.5 

287.3 

267.3 

151.3 

295.3 

347.1 

Sensitivity  tests  with  constant  bottom  friction  coefficients  of  0.002,  0.003,  and 
0.004  showed  that  the  differences  in  the  tidal  elevation  and  the  phase  lag  were  hardly 
noticeable  (case  6,  case  7,  case  8 in  Table  5.5,  Table  5.6,  Table  5.7).  Phase  lag 
increased  slightly  with  the  increased  bottom  friction  coefficient.  When  the  law  of 
the  wall  was  used  (case  3),  the  tidal  elevation  and  phase  lag  were  improved  with 
respect  to  the  observation  data  when  compared  to  results  from  those  obtained  with 
constant  bottom  friction  coefficients.  However,  the  tidal  elevation  at  Station  1 was 
overestimated  both  with  the  constant  bottom  friction  coefficients  and  with  the  law  of 
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the  wall.  When  the  Monin- Oboukhov  length  was  used  with  the  law  of  the  wall,  the 
result  was  improved  substantially  (case  2).  Sensitivity  to  bottom  roughness  height 
was  performed  with  0.4,  0.1,  and  0.04  in  case  2,  case  4,  and  case  5,  respectively. 
With  the  smaller  roughness  height,  the  tidal  elevation  increased  and  the  phase  lag 
decreased.  In  the  fine  grid  system,  case  9 and  case  10  showed  the  same  trend  for 
the  tidal  elevation  as  in  coarse  grid  system  for  the  test  of  roughness  height.  For  u- 
velocity,  u- velocity,  both  amplitude  and  phase  lag  increased  with  a small  roughness 
height  (case  9 and  case  10  in  Table  5.8  and  Table  5.9). 

5.3.2  Effect  of  Bathymetry 

Sensitivity  of  model  prediction  to  the  bathymetry  was  examined  with  several 
different  sets  of  bathymetry  used  in  the  simulations.  Naturally  bathymetry  is  one 
of  the  most  important  simulation  parameters  responsible  for  accurate  model  results. 
The  bathymetry  change  is  directly  related  to  the  bottom  friction.  Depth  data  were 
corrected  many  times  to  best  represent  the  depth  near  each  grid  point.  Correction 
of  depth  data  in  a certain  grid  system  must  be  performed  in  the  model  calibration 
process  in  order  to  use  the  model  as  prediction  tool.  Also  experience  shows  that  a 
model  that  does  not  accurately  resolve  the  deeper  channel  of  an  estuary  system  cannot 
predict  the  correct  degree  of  stratification  and  two-layer  estuarine  residual  flow. 

5.3.3  Effect  of  Horizontal  Grid 

The  discretization  of  the  governing  equations  leads  unavoidably  to  approximation 
errors.  Two  grid  systems:  one  coarse  grid  (Figure  5.1)  and  one  fine  grid  (Figure  5.2) 
were  used  in  this  study.  Compared  to  the  coarse  grid  system,  the  grid  number  in  the 
transverse  direction  is  doubled  in  fine  grid  system:  coarse  grid  system  has  only  4 to 
5 grid  lines  in  tj  direction,  while  fine  grid  system  has  8 to  10  grid  lines.  Case  2 and 
case  10  show  the  deviations  in  elevations  and  velocities  between  two  grid  systems. 
As  shown  in  Tables  5.5,  5.6,  and  5.7,  it  is  difficult  to  discern  the  differences  between 
the  coarse-grid  and  fine-grid  results  at  stations  2,  3,  and  4.  However  in  the  fine  grid 
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Table  5.7:  Comparison  of  data  and  model  tidal  con- 
stituents of  surface  u-velocity  at  Station  4 (amplitude 
in  cm,  phase  lag  in  degree,  6/19  - 7/17/85). 


constituent 

M2 

S2 

N2 

K1 

01 

data 

amplitude 

51.5 

8.8 

10.0 

7.1 

5.5 

phase  lag 

122.72 

137.76 

355.16 

110.50 

146.66 

case9 

amplitude 

52.4 

7.7 

8.3 

3.3 

4.0 

phase  lag 

102.94 

81.50 

321.81 

63.29 

130.20 

case 10 

amplitude 

53.5 

8.7 

9.1 

4.6 

5.2 

phase  lag 

107.48 

84.76 

328.88 

66.85 

131.18 

easel  1 

amplitude 

52.1 

8.8 

9.2 

4.7 

5.7 

phase  lag 

109.77 

88.73 

333.35 

69.54 

127.95 

case 12 

amplitude 

52.2 

7.5 

8.3 

3.6 

4.1 

phase  lag 

102.88 

80.24 

322.65 

64.81 

126.44 

system,  the  circulation  gyre  near  the  Hampton  Roads  (Figure  5.26  and  Figure  5.27) 
is  much  clearer.  The  two  grid  systems  represent  very  similar  physical  processes. 

5.3.4  Effect  of  Time  Step 

Various  time  steps  of  120,  180  and  360  seconds  were  used  with  a bottom  friction 
factor  0.003  in  the  coarse  grid  system.  The  deviations  of  the  elevations  and  the 
velocities  for  each  selected  time  step  at  the  selected  stations  were  negligible.  The 
results  show  that  time  step  interval  is  not  a restriction  for  stability  because  this 
model  was  run  with  implicit  method. 

5.4  Neap-Spring  Stratification  and  Destratification  Cycle 

In  a partially-mixed  estuary  like  James  River,  the  degree  of  stratification  and 
the  circulation  pattern  are  controlled  by  the  volume  of  water  flowing  up  the  estuary 
during  the  flood  tide  and  the  volume  of  fresh  water  flowing  down  the  estuary.  A 
primary  source  of  mixing  energy  is  tidal  forcing.  The  effect  associated  with  neap-to- 
spring  tide  in  James  River  is  seen  in  Figure  5.12  and  Figure  5.13  near  the  James  River 
bridge.  During  the  spring  tide  with  high  tide  range  (July  3),  well  mixed  condition 
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Table  5.8:  Comparison  of  data  and  model  tidal  con- 
stituents of  surface  v-velocity  at  Station  4 (amplitude 
in  cm,  phase  lag  in  degree,  6/19  - 7/17/85). 


constituent 

M2 

S2 

N2 

K1 

01 

data 

amplitude 

39.5 

2.7 

5.8 

5.5 

3.8 

phase  lag 

302.34 

324.92 

176.36 

296.27 

321.37 

case9 

amplitude 

36.6 

4.5 

4.4 

1.9 

2.8 

phase  lag 

298.04 

275.47 

152.24 

274.85 

326.79 

case 10 

amplitude 

38.3 

5.3 

5.4 

2.6 

3.8 

phase  lag 

299.29 

277.20 

158.76 

275.72 

323.91 

easel 1 

amplitude 

37.7 

5.4 

5.6 

2.6 

4.0 

phase  lag 

299.98 

280.69 

160.17 

272.97 

320.12 

easel 2 

amplitude 

36.4 

4.6 

4.5 

2.0 

2.8 

phase  lag 

297.84 

272.06 

151.59 

275.04 

325.32 

was  developed,  while  during  the  ebb  tide  with  lower  tide  range  (July  9),  stratification 
became  more  apparent. 
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Figure  5.12:  Simulated  longitudinal- vertical  salinity  at  slackwater  before  flood  on 
July  3,  1985. 
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Figure  5.13:  Simulated  longitudinal-vertical  salinity  at  slackwater  before  flood  on 
July  9,  1985. 
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The  degree  of  stratification  can  be  related  to  buoyancy  and  mixing  energy  through 
the  non-dimensional  parameter  M (Hansen  and  Rattray,  1965)  which  was  defined  in 
section  3.5.  As  M increases,  stratification  decreases  until  it  becomes  well-mixed. 
However,  it  is  difficult  to  evaluate  the  mixing  parameters  and  relate  them  to  tidal 
action.  Cerco  (1982)  used  an  alternative  dimensionless  quantity  which  expresses  the 
ratio  of  buoyancy  to  mixing  action  with  uj/ut,  i.e.,  the  ratio  of  freshwater  velocity  to 
tidal  velocity.  He  divided  two  regions  delineated  by  A s/s  = 0.1  and  uj/ut  = 0.008. 
For  values  of  Uf/ut  > 0.008,  observed  stratification  is  generally  greater  than  10%.  For 
values  of  Uf/ut  < 0.008,  stratification  is  generally  less  than  10%.  Spring  tides  produce 
stronger  tidal  currents,  and  hence  more  intense  mixing  than  neap  tides.  According  to 
his  analysis  in  the  James  River  data,  the  cyclical  nature  of  stratification  in  the  estuary 
and  the  relationship  of  the  stratification-destratification  cycle  to  the  spring-neap  tidal 
cycle  are  clear. 

Model  response  to  variation  in  tidal  range  was  examined  through  the  29-day 
simulation  of  the  baseline  run  (case  2).  Figures  5.14,  5.15,  and  5.16  show  the  salinity 
and  velocity  time-series  during  the  29-day  period  at  three  levels  (surface,  mid,  and 
bottom  level)  at  Station  1,  Station  2,  and  Station  3,  respectively.  The  effect  of  the 
spring-neap  tidal  cycle  was  apparent  in  the  time  series  of  salinity  and  velocity.  An 
increase  in  tidal  range  during  the  spring  tide  results  in  an  increase  in  velocity  and 
a decrease  in  stratification  while  a decrease  in  tidal  range  results  in  an  increase  in 
stratification. 

The  near-surface  longitudinal  salinity  distributions  on  July  3 (spring  tide),  on 
July  9 (neap  tide),  and  July  17  (spring  tide)  during  the  29-day  simulation  are  shown 
in  Figures  5.17,  5.18,  and  5.19,  respectively.  The  measured  profile  is  generally  within 
the  range  of  the  simulated  salinity  profiles.  The  range  of  variation  in  time  changes 
according  to  the  spring-neap  tidal  variation.  It  is  clear  that  the  vertical  salinity 
profiles  go  through  cycles  of  stratification  and  destratification  in  direct  response  to 
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Figure  5.14:  Computed  elevation,  and  velocity  and  salinity  at  surface,  mid,  and 
bottom  level  for  29-day  simulation  at  Station  1. 


146 


JUNE  1 9 - JULY  17, 1985  AT  (65,1 8) 


500  600 

Time  (Hours) 


Figure  5.15:  Computed  elevation,  and  velocity  and  salinity  at  surface,  mid,  and 
bottom  level  for  29-day  simulation  at  Station  2. 
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Figure  5.16:  Computed  elevation,  and  velocity  and  salinity  at  surface,  mid,  and 
bottom  level  for  29-day  simulation  at  Station  3. 
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Longitudinal  salinity  profile  on  July  3, 1985 


Figure  5.17:  Model  predicted  longitudinal  surface  salinity  distribution  with  3000 
ft3 /sec  of  river  discharge  on  July  3,  1985. 
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Figure  5.18:  Model  predicted  longitudinal  surface  salinity  distribution  with  3000 
ft3 /sec  of  river  discharge  on  July  9,  1985. 
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Figure  5.19:  Model  predicted  longitudinal  surface  salinity  distribution  with  3000 
ft3 /sec  of  river  discharge  on  July  17,  1985. 
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the  varying  tidal  conditions  and  vertical  mixing  in  the  water  column. 

Model  results  corresponding  to  the  slackwater  salinity  distributions  on  July  3 and 
July  9,  1985  are  presented  in  Figures  5.12  and  5.13.  These  model  results  exhibit  a sig- 
nificant spring- to-neap  variation  in  salinity  distribution  as  was  found  in  the  measured 
data. 

5.5  Salinity  Intrusion 

The  salinity  intrusion  in  an  estuary  is  subjected  to  longitudinal  displacement 
and  changes  in  vertical  structure  as  a result  of  variations  in  freshwater  runoff  and 
tidal  motion.  The  location  of  the  upstream  limit  of  the  salinity  intrusion  is  deter- 
mined primarily  by  the  force  balance  between  the  downstream  hydrostatic  head  of  the 
freshwater  flow  and  the  upstream  pressure  gradient  due  to  the  longitudinal  salinity 
gradient.  The  intrusion  moves  downstream  with  increasing  freshwater  runoff  and  up- 
stream with  decreasing  runoff.  In  the  James  River,  during  summer  and  fall,  low-flow 
conditions  prevail  and  the  intrusion  gradually  moves  upstream  until  the  high  runoff 
during  winter  and  spring  recurs  (Cerco,  1982).  In  this  study,  because  of  the  lack  of 
field  data,  model  validation  using  time- varying  river  discharge  was  not  made.  Instead, 
a series  of  sensitivity  tests  with  constant  river  discharge  was  performed  to  investigate 
the  simulated  response  of  the  estuary  to  river  flow.  Figures  5.17-5.25  show  the  lon- 
gitudinal salinity  distribution  at  the  surface  level  (K=ll)  along  the  channel  (K=18) 
on  July  3 (spring),  on  July  9 (neap),  and  on  July  17  (spring),  with  a constant  river 
flow  rate  of  3000  ft3/s.  Each  figure  contains  6 salinity  distribution  lines  plotted  with 
two  hour  increaments.  The  shapes  of  the  curves  are  approximately  described  by  the 
the  hyperbolic  tangent  function  which  is  common  for  many  coastal-plain  estuaries. 

Hansen  and  Rattray  (1965)  developed  a pair  of  steady-state  stream-function  equa- 
tions describing  circulation  and  salt  balance  in  an  estuary.  They  solved  the  equations 
by  defining  three  salinity  regimes:  outer  regime,  central  regime,  and  inner  regime.  In 
the  central  regime,  the  vertical  salinity  stratification  is  nearly  independent  of  longi- 
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Figure  5.20:  Model  predicted  longitudinal  surface  salinity  distribution  with  1000 
ft3 /sec  of  river  discharge  on  July  3,  1985. 
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Figure  5.21:  Model  predicted  longitudinal  surface  salinity  distribution  with  1000 
ft3/sec  of  river  discharge  on  July  9,  1985. 
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Figure  5.22:  Model  predicted  longitudinal  surface  salinity  distribution  with  1000 
ft3 /sec  of  river  discharge  on  July  17,  1985. 
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Figure  5.23:  Model  predicted  longitudinal  surface  salinity  distribution  with  9000 
ft3 /sec  of  river  discharge  on  July  3,  1985. 
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Figure  5.24:  Model  predicted  longitudinal  surface  salinity  distribution  with  9000 
ft3/sec  of  river  discharge  on  July  9,  1985. 
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Figure  5.25:  Model  predicted  longitudinal  surface  salinity  distribution  with  9000 
ft3/ sec  of  river  discharge  on  July  17,  1985. 
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tudinal  position.  In  the  inner  and  outer  regimes,  the  degree  of  stratification  becomes 
smaller  in  proportion  to  the  distance  from  the  central  regime.  However  the  shape  of 
the  longitudinal  salinity  profile  constantly  moves  up  and  down  the  estuary  in  response 
to  changes  in  river  discharge.  Figures  5.17-5.25  also  show  that  the  vertical  range  of 
salinity  variation  changes  in  accordance  with  the  spring-neap  variation.  As  expected, 
during  spring  tide  the  variation  of  salinity  distribution  is  large,  and  the  thickness, 
hence,  is  large. 

Three  additional  river  discharges  were  tested:  1000  ft3/sec  (Figures  5.20-5.22) 
which  corresponds  to  low  flow  conditions,  3000  ft3/sec  (Figures  5.17-5.19)  which 
corresponds  to  the  average  flow  condition  during  the  29-day  simulation  starting  June 
19,  1985,  and  9000  ft3 /sec  (Figures  5.23-5.25)  which  is  a usual  mean  annual  value. 
With  the  largest  river  discharge  of  9000  ft3 /sec,  the  head  of  salinity  intrusion  was 
found  near  70  km  from  the  mouth.  With  the  smallest  value  of  river  discharge  of  1000 
ft3 /sec,  the  head  of  salinity  intrusion  was  found  near  90  km  upstream.  With  3000 
ft3 /sec  of  river  discharge,  the  head  of  salt  intrusion  was  located  around  80  km  from 
the  river  mouth. 

In  the  model  simulation,  concurrent  with  the  high  flow  rate,  the  intrusion  is 
pushed  to  its  most  downstream  location.  With  low  flow  rate,  the  intrusion  moves 
upstream.  The  salinity  variation  with  time  is  larger  at  spring  tide  than  at  neap  tide. 
Higher  runoff  also  causes  a larger  salinity  variation  with  time. 

5.6  Frontal  Dynamics 

To  examine  this  frontal  structure  in  the  model  simulation,  the  model  was  run 
for  29  days  starting  from  June  19,  1985  in  the  fine  grid  system  (case  10).  In  the 
coarse  grid  system,  it  was  hard  to  see  a gyre  or  a convergence  area  because  only  4 
or  5 grid  cells  were  constructed  in  the  transverse  direction  and  it  was  impossible  to 
resolve  the  real  depth  configuration.  In  the  fine  grid  system,  the  number  of  cells  in  the 
transverse  direction  was  8 or  9.  With  these  grid  cells,  we  could  see  a gyre  because  the 
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deep  channel  was  resolved.  The  model  results  on  July  3,  1985  were  taken  and  plotted 
for  one  tidal  cycle.  Time  step  used  here  was  360  seconds.  Hourly-plottings  were  made 
for  instantaneous  velocity  field  and  salinity  contour  at  the  surface  (Figures  5.26-5.36), 
and  for  instantaneous  velocity  flield  and  salinity  contour  at  bottom  (Figures  5.37- 
5.47).  Velocity  structures  in  u and  w direction  are  also  plotted  along  the  channel 
(Figures  5.48-5.58)  near  Newport  News  and  along  the  one-grid  away  from  the  channel 
(Figures  5.59-5.69),  and  time  evolution  of  combined  u,  w,  salinity,  and  elevation  along 
the  one-grid  away  from  the  channel  (Figures  5.70-5.80),  where  u is  at  surface,  w is  at 
mid-depth,  and  salinity  is  at  surface,  mid,  and  bottom. 

Two  factors  are  essential  to  the  formation  of  the  front.  They  are  the  convergence 
of  surface  currents  and  the  density  contrast  between  the  convergent  flows.  As  the 
heavier  water  pushes  against  the  lighter  water  at  the  convergence,  it  dives  beneath 
the  lighter  water  as  a result  of  gravitational  force.  On  the  flood  phase  of  tide,  the 
incoming  water  in  Hampton  Roads  meets  the  ebbing  water  of  the  lower  James  off 
Newport  News  Point  (Figures  5.26-5.28).  Especially  a gyre  was  formed  at  hour  1 
(Figures  5.26).  Bottom  currents  also  show  a gyre  circulation  (Figures  5.37).  This 
convergence  zone  near  Newport  News  Point  was  maintained  about  3 hours  after  the 
flood  tide  started  at  the  mouth.  The  density  difference  between  the  two  water  masses 
is  large  enough  to  make  the  more  saline  Hampton  Roads  water  dive  beneath  the  fresh 
water  of  the  lower  James.  When  the  salinity  contour  plots  are  compared  with  velocity 
fields  (Figures  5.26-5.36  and  Figures  5.37-5.47),  the  salinity  contour  lines  move  up 
and  down  according  to  velocity  field.  Therefore,  the  advection  of  the  salinity  plays  an 
apparently  important  role  for  the  salinity  distribution.  As  the  heavier  water,  either 
from  Hampton  Flats/Newport  News  Bar  or  flowing  along  Newport  News  Channel, 
meets  the  lighter  water  off  Newport  News  Point,  it  also  encounters  an  abrupt  change 
of  bottom  topography.  This  steep  topographic  gradient  not  only  enhances  the  diving 
of  heavier  water  but  also  tends  to  stabilize  the  location  of  the  front.  Figures  5.59-5.62 
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show  a clear  picture  of  the  topographical  effects.  We  can  see  a plume  type  of  current 
movement  near  the  abrupt  change  of  bottom  topography.  Also  we  can  see  a strong 
vertical  transport.  For  the  period  of  frontal  zone  formation,  a gyre  of  vertical  eddy 
motion  appeared  because  of  bottom  topography  and  strong  convergence  of  ebb  and 
flood  (Figures  5.50  and  Figures  5.61).  Figures  5.70-5.80  show  the  variation  of  salinity 
structure  and  vertical  velocity  during  the  front  formation.  The  maximum  downward 
vertical  velocity  was  seen  at  hour  3 (Figure  5.72). 

In  the  model  simulation,  strong  vertical  velocities  are  shown  near  the  SBF  through 
the  front  moving  upestuary  with  time  (Figures  5.59-5.62). 

Results  of  current  measurements  indicate  that  flood  tide  currents  on  Hampton 
Flats  lead  the  current  along  the  Newport  News  Channel  section.  This  is  shown  in 
the  model  simulation.  Figure  5.26  shows  a recirculation  eddy  on  the  east  end  of 
Hampton  Flats  while  the  current  in  the  channel  is  still  ebbing.  This  phenomenon 
was  simulated  in  the  hydraulic  model  test  by  Byrne  (1987).  When  the  currents 
within  the  channel  section  are  slack,  the  flood  current  over  the  shallow  Hampton 
Flats  converges  at  Newport  News  Point  (Figure  5.26).  There  it  encounters  ebbing 
water  of  less  density  than  the  lower  and  dives  beneath  the  less  dense  water.  As  the 
front  moves  over  the  topographic  gradient  at  the  western  end  of  the  dredged  channel, 
the  increase  in  water  depth  slows  down  the  upriver  movement  of  the  front  and  also 
increases  the  diving  action  of  the  surface  water  (Figure  5.27).  The  substantial  lead 
of  flood  current  on  Hampton  Flats  over  those  in  deep  waters  in  the  channel  may  be 
attributed  to  several  mechanisms.  The  geometry  near  Newport  News  Point  develops 
a topographic  eddy,  where  the  channel  makes  a 90°  bend.  As  the  ebb  flow  from  the 
lower  James  tries  to  bend,  a topographic  eddy  begins  to  evolve  on  the  shallow  flanks. 
Toward  the  late  stage  of  ebb  tide,  the  eddy  is  fully  developed  into  a counterclockwise 
circulation  over  Hampton  Flats.  The  circulation  advances  the  flood  phase  of  tide 
and  enhances  the  strength  of  flood  currents  in  its  early  stage.  One  dimensional 
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analysis  of  tidal  dynamics  indicates  that  the  friction  tends  to  advance  the  phase  of 
tidal  current  relative  to  that  of  tidal  stage  (Ippen,  1966,  p.509).  Boundary  layer 
theory  also  indicates  that  the  oscillating  current  near  the  boundary  will  lead  that 
of  the  interior  flow  (Lamb,  1962,  p.620).  Therefore  it  is  expected  that  flood  current 
on  Hampton  Flats  would  lead  those  in  the  channel.  The  friction  due  to  both  the 
shallow  depth  and  lateral  boundary,  plays  a role  in  the  phase  variation  of  tidal  current 
around  Newport  News  Point.  Another  mechanism  in  phase  sequence  is  a damped  co- 
oscillating tide.  Newport  News  Point  is  about  one-half  tidal  wave  length  from  the 
head  of  tide.  Therefore  anti-node  of  a damped  co-oscillating  tide  is  located  at  this 
location.  Tide  table  data  show  a local  maximum  of  tidal  range  at  this  location,  and  a 
local  minimum  at  an  upriver  location  about  one-quarter  wave  length  from  the  head  of 
tide.  A damped  co-  oscillating  tide  would  possess  the  same  features.  The  flood  current 
on  the  downriver  side  starts  much  earlier  than  that  on  the  upriver  side,  resulting  in 
convergent  flow  at  the  anti-node  during  the  early  stage  of  flood  tide. 

The  density  field  exhibits  the  usual  trends  of  a coastal  plain  estuary,  with  moder- 
ate density  increases  with  depth  and  decreases  with  distance  from  the  mouth.  Obser- 
vations have  revealed  that  these  trends  are  persistent  features  with  highly  localized, 
sharp  horizontal  gradients  (fronts)  in  the  near-surface  density,  occurring  systemati- 
cally at  particular  locations  and  during  a particular  phase  of  the  tidal  cycle.  Early 
in  the  flood,  a region  of  relatively  strong  vertical  salinity  gradient  developed  in  the 
upriver  side  front.  Also  a region  of  strong  horizontal  gradient  is  found  near  the  front. 
With  time  evolution,  the  front  has  moved  upriver.  These  density  field  variations  in 
the  model  simulation  agree  with  the  observed  data  by  Byrne  (1987)  (Figures  5.26- 
5.36).  In  the  longitudinal  sections,  in  addition  to  the  horizontal  transition  of  the 
surface  front,  there  is  also  an  apparent  deepening  and  steepening  of  the  frontal  zone 
as  the  flood  tide  progresses  (Figures  5.81-5.91).  However,  this  is  not  a straightforward 
change  of  position  and  shape  of  a well-defined  interface.  Rather,  higher  salinity  water 
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has  been  incorporated  into  the  downriver  side  of  the  high-gradient  region,  but  the 
depth  at  which  each  isohaline  approaches  the  horizontal  on  the  upriver  side  remains 
nearly  constant  throughout  the  tidal  cycle. 
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Figure  5.26:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  1 a.m.  1 on  July  3,  1985. 
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Figure  5.27:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  2 a.m.  on  July  3,  1985. 
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Figure  5.28:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  3 a.m.  on  July  3,  1985. 
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Figure  5.29:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  4 a.m.  on  July  3,  1985. 
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Figure  5.30:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  5 a.m.  on  July  3,  1985. 
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Figure  5.31:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  6 a.m.  on  July  3,  1985. 
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Figure  5.32:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  7 a.m.  on  July  3,  1985. 
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Figure  5.33:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  8 a.m.  on  July  3,  1985. 
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Figure  5.34:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  9 a.m.  on  July  3,  1985. 
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Figure  5.35:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  10  a.m.  on  July  3,  1985. 
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Figure  5.36:  Instantaneous  near-surface  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  11  a.m.  on  July  3,  1985. 
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Figure  5.37:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  1 a.m.  on  July  3,  1985. 
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Figure  5.38:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  2 a.m.  on  July  3,  1985. 
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Figure  5.39:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  3 a.m.  on  July  3,  1985. 
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Figure  5.40:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  4 a.m.  on  July  3,  1985. 
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Figure  5.41:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  5 a.m.  on  July  3,  1985. 
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Figure  5.42:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  6 a.m.  on  July  3,  1985. 
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Figure  5.43:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  7 a.m.  on  July  3,  1985. 
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Figure  5.44:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  8 a.m.  on  July  3,  1985. 
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Figure  5.45:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  9 a.m.  on  July  3,  1985. 
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VARIABLE:  velocity  and  salinity 
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Figure  5.46:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  10  a.m.  on  July  3,  1985. 
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Figure  5.47:  Instantaneous  near-bottom  horizontal  velocity  and  salinity  ( ppt ) distri- 
bution in  lower  James  River  at  11  a.m.  on  July  3,  1985. 
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Figure  5.48:  Instantaneous  velocity  in  ( x,z ) domain  along  the  channel  in  lower  James 
River  at  1 a.m.  on  July  3,  1985. 
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Figure  5.49:  Instantaneous  velocity  in  (x,z)  domain  along  the  channel  in  lower  James 
River  at  2 a.m.  on  July  3,  1985. 
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Figure  5.50:  Instantaneous  velocity  in  (x,  z)  domain  along  the  channel  in  lower  James 
River  at  3 a.m.  on  July  3,  1985. 
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Figure  5.51:  Instantaneous  velocity  in  (x,z)  domain  along  the  channel  in  lower  James 
River  at  4 a.m.  on  July  3,  1985. 


189 


Figure  5.52:  Instantaneous  velocity  in  (x,z)  domain  along  the  channel  in  lower  James 
River  at  5 a.m.  on  July  3,  1985. 
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Figure  5.53:  Instantaneous  velocity  in  (x,  z ) domain  along  the  channel  in  lower  James 
River  at  6 a.m.  on  July  3,  1985. 
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Figure  5.54:  Instantaneous  velocity  in  (x,z)  domain  along  the  channel  in  lower  James 
River  at  7 a.m.  on  July  3,  1985. 
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Figure  5.55:  Instantaneous  velocity  in  ( x,z ) domain  along  the  channel  in  lower  James 
River  at  8 a.m.  on  July  3,  1985. 
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Figure  5.56:  Instantaneous  velocity  in  (x,z)  domain  along  the  channel  in  lower  James 
River  at  9 a.m.  on  July  3,  1985. 
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Figure  5.57:  Instantaneous  velocity  in  (x,  z)  domain  along  the  channel  in  lower  James 
River  at  10  a.m.  on  July  3,  1985. 


195 


Figure  5.58:  Instantaneous  velocity  in  ( x , z)  domain  along  the  channel  in  lower  James 
River  at  11  a.m.  on  July  3,  1985. 
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Figure  5.59:  Instantaneous  velocity  in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  1 a.m.  on  July  3,  1985. 
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Figure  5.60:  Instantaneous  velocity  in  (x,z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  2 a.m.  on  July  3,  1985. 
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Figure  5.61:  Instantaneous  velocity  in  ( x,z ) domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  3 a.m.  on  July  3,  1985. 
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Figure  5.62:  Instantaneous  velocity  in  ( x,z ) domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  4 a.m.  on  July  3,  1985. 
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Figure  5.63:  Instantaneous  velocity  in  (x,z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  5 a.m.  on  July  3,  1985. 
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Figure  5.64:  Instantaneous  velocity  in  (x,z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  6 a.m.  on  July  3,  1985. 
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Figure  5.65:  Instantaneous  velocity  in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  7 a.m.  on  July  3,  1985. 


203 


Figure  5.66:  Instantaneous  velocity  in  ( x,z ) domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  8 a.m.  on  July  3,  1985. 
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Figure  5.67:  Instantaneous  velocity  in  ( x,z ) domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  9 a.m.  on  July  3,  1985. 
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Figure  5.68:  Instantaneous  velocity  in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  10  a.m.  on  July  3,  1985. 
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Figure  5.69:  Instantaneous  velocity  in  ( x,z ) domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  11  a.m.  on  July  3,  1985. 
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SALINITY  AND  VELOCITY  NEAR  THE  FRONT 


Figure  5.70:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  1 a.m.  on  July  3,  1985. 
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SALINITY  AND  VELOCITY  NEAR  THE  FRONT 


Figure  5.71:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  2 a.m.  on  July  3,  1985. 
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SALINITY  AND  VELOCITY  NEAR  THE  FRONT 


Figure  5.72:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  3 a.m.  on  July  3,  1985. 
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Figure  5.73:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  4 a.m.  on  July  3,  1985. 
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SALINITY  AND  VELOCITY  NEAR  THE  FRONT 


TIME (HOUR)  : 5 


Figure  5.74:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  5 a.m.  on  July  3,  1985. 
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Figure  5.75:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  6 a.m.  on  July  3,  1985. 


W-VELOC I TY( CM/SEC) 


SURFACE  ELEVATION (CM) 


213 


SALINITY  AND  VELOCITY  NEAR  THE  FRONT 


Figure  5.76:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  7 a.m.  on  July  3,  1985. 
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SALINITY  AND  VELOCITY  NEAR  THE  FRONT 


Figure  5.77:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  8 a.m.  on  July  3,  1985. 
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SALINITY  AND  VELOCITY  NEAR  THE  FRONT 


Figure  5.78:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  9 a.m.  on  July  3,  1985. 
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Figure  5.79:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  10  a.m.  on  July  3,  1985. 
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SALINITY  AND  VELOCITY  NEAR  THE  FRONT 


Figure  5.80:  u at  surface,  w at  mid-depth,  salinity  at  surface-,  mid-,  and  bot- 
tom-depth, and  elevation  along  the  one-grid  away  from  the  channel  in  lower  James 
River  at  11  a.m.  on  July  3,  1985. 
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Figure  5.81:  Salinity  contours  ( ppt ) in  ( x,z ) domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  1 a.m.  on  July  3,  1985. 
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Figure  5.82:  Salinity  contours  ( ppt ) in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  2 a.m.  on  July  3,  1985. 
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Figure  5.83:  Salinity  contours  ( ppt ) in  (x,z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  3 a.m.  on  July  3,  1985. 
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Figure  5.84:  Salinity  contours  ( ppt ) in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  4 a.m.  on  July  3,  1985. 
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Figure  5.85:  Salinity  contours  ( ppt ) in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  5 a.m.  on  July  3,  1985. 
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Figure  5.86:  Salinity  contours  ( ppt ) in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  6 a.m.  on  July  3,  1985. 
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Figure  5.87:  Salinity  contours  ( ppt ) in  (x,z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  7 a.m.  on  July  3,  1985. 
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Figure  5.88:  Salinity  contours  ( ppt ) in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  8 a.m.  on  July  3,  1985. 
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Figure  5.89:  Salinity  contours  ( ppt ) in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  9 a.m.  on  July  3,  1985. 
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Figure  5.90:  Salinity  contours  ( ppt ) in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  10  a.m.  on  July  3,  1985. 
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Figure  5.91:  Salinity  contours  ( ppt ) in  (x,  z)  domain  along  the  one-grid  away  from 
the  channel  in  lower  James  River  at  11  a.m.  on  July  3,  1985. 
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5.7  Model  Calibration  and  Validation 

Usually  model  performance  has  been  assessed  based  on  the  subjective  inspection 
of  plotting  results.  The  objective  analysis  of  model  performance  is  important  because 
it  tells  model  skill  quantitatively  and  makes  it  possible  to  compare  one  model  to  oth- 
ers. NOS  recently  proposed  the  standards  for  tides  and  currents  as  an  assessment  of 
model  accuracy.  To  quantify  differences  between  model  output  and  observed  data, 
Hess  (1991)  performed  statistical  analysis  of  model  values  and  observed  values,  in- 
cluding rms  differences,  gain,  mean  time  lag,  and  rms  time  lag.  Sheng  and  Peene 
(1992)  used  similar  statistical  analysis  to  compare  their  model  results  with  observed 
water  level  and  current  in  Tampa  Bay  and  Sarasota  Bay. 

First  consider  a series  of  demeaned  observed  values,  y,  (for  i = 1,  • • • , n),  and  a 
series  of  demeaned  model  output  values,  y,-,  with  the  same  starting  time,  the  same 
number  of  values,  and  uniform  time  intervals.  The  difference  between  values  in  the 
modeled  and  observed  time  series  is 

di  = y,  - Vi  (5.5) 

and  the  rms  difference  is 

Drms  = EK-)7"]1/2  (5.6) 

i=i 

Since  tidal  forcing  is  important  in  most  estuaries,  statistics  for  tidal  extrema  (highs 
and  lows,  floods  and  ebbs)  are  used  for  the  statistical  analysis.  The  observed  extrema 
Yj  at  times  7),  for  j = 1,2 , • • • , TV  and  the  model  extrema  Yj  at  times  Tj,  for  j — 
1, 2,  • • • , N are  found.  The  mean  tidal  range,  R is 

= (5-7) 

j=  1 

Range  is  defined  as  difference  from  maximum  to  minimum.  Since  observed  tide  ranges 
can  vary  widely  from  station  to  station,  a more  meaningful  statistic  is  the  ratio  of  the 
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rms  difference  to  the  mean  range  (/?), 

D'  = Drms/R  (5.8) 

The  observed  and  modeled  extrema  are  then  compared  by  selecting  for  each  individual 
observed  extremum  and  the  closest  modeled  extremum.  There  will  be  v pairs.  A gain 
( G ),  ratio  of  modeled  to  observed  extrema  amplitude,  is  defined  as 

o = : ifiriYj >/*  (5.9) 

J=1 

When  the  observed  tide  range  is  non-uniform  and  small,  a weighted  gain  ( Gw ) is  used 
with  the  weight  being  equal  to  the  individual  modeled  range  divided  by  R. 

G»  = (l)'t\Yi\/N  (5.10) 

n 3=1 

The  weighted  gain  describes  characteristics  of  the  whole  series.  The  rms  differences 
for  extrema  (Erms)  gives  the  information  on  individual  differences,  which  can  be 
expressed, 

Er*.  = E(*5>  - YtfW*  (5.11) 

3= 1 

For  phase  differences,  a mean  time  lag  (Lm)  for  the  extrema  is  defined  as 

Lm  = t(T,  - Tj)/v  (5.12) 

i=i 

The  rms  lag  (Lrma)  is  defined  as 


in».  = E (Ti'-Tj)2/*]'!1  (5.13) 

3= 1 

The  Lm  can  show  the  modeled  wave  speed  errors  as  a lag  (or  lead)  at  all  stations.  But 
the  Lrms  is  valuable  when  Lm  is  small.  When  the  results  from  stations  are  compared, 
these  multistation  summary  statistics  are  referred  to  as  global  quantities, 

i= 1 i=l 


(5.14) 
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where  v is  the  number  of  data  pairs,  h is  the  sampling  interval  for  each  station,  and 
s is  the  number  of  stations. 

Now  the  ratio  of  the  Erms  difference  to  the  mean  range  ( R ) and  the  ratio  of  LTTns 
to  6.21  hours  (one  fourth  of  M2  tidal  period)  are  used  to  calculate  the  model  skill 
parameters  ( E'  = Erms/R  and  L'  = Lrms/ 6.21).  L'  approximates  the  maximum 
time  between  extrema.  Three  skill  parameters  are  obtained,  which  are  Sp  = 1 — D' , 
Se  — 1 — E',  and  Si  — 1 — L' . These  skill  parameters  can  be  used  to  judge  model 
accuracy. 

During  the  sensitivity  analysis,  model  inputs  are  varied  to  assess  changes  in  the 
solution  for  the  parameters  of  interest.  The  results  of  a 29-day  run  for  water  levels 
and  currents  are  shown  in  Table  5.9,  Table  5.10,  Table  5.11,  and  Table  5.12.  Case 
2 was  chosen  as  a calibrated  result,  which  used  the  law  of  the  wall  with  the  Monin- 
Oboukhov  length  with  bottom  roughness  height  0.4  cm.  Hourly  rms  amplitude  dif- 
ferences (Z)rma=0.03)  are  small,  model  gain  ( Gw  = 1.01)  is  nearly  unity,  and  rms  time 
differences  ( Lrms  = 30 minutes)  are  relatively  high.  Mean  lag  ( Lm  = —4 minutes)  is 
small  and  suggests  no  bias  in  propagation  speed  errors.  Model  skill  parameters  are 
all  greater  than  95  %. 

Statistics  for  water  currents  from  the  same  run  of  case  2 are  shown  in  Ta- 
ble 5.11.  Model  gain  for  currents  ( Gw  = 0.87)  showed  that  the  currents  were  un- 
derestimated. The  rms  time  differences  ( LTms  = 36  minutes)  was  large,  and  mean 
lag  ( Lm  = 22  minutes)  was  worse  compared  to  that  in  the  water  elevation.  Model 
skill  parameters  are  about  90  % (SD  = 88%,  SE  = 90%,  SL  = 90%).  The  model 
can  be  improved  by  the  increase  of  the  model  grid  resolution  and  the  more  accurate 
bottom  topography. 

The  results  from  the  calibration  runs  show  that  all  the  individual  model  skill  parame- 
ters for  the  surface  elevations  are  greater  than  90  %,  so  small  changes  in  the  roughness 
height  and  bottom  friction  coefficient  are  not  so  significant. 
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Table  5.9:  Model  skill  parameters  for  the  sensitivity  test  to  surface  elevation  at  Station 
1.  


R 

Drms 

Drms/  R 

V 

Gw 

E 

J-/ rms 

Lm 

Lrms 

easel 

65.85 

4.78 

0.07 

111 

1.18 

6.59 

-0.03 

0.37 

case2 

65.85 

2.46 

0.04 

111 

1.04 

2.27 

-0.09 

0.36 

case3 

65.85 

4.19 

0.06 

111 

1.17 

5.59 

-0.03 

0.37 

case4 

65.85 

2.64 

0.04 

111 

1.09 

3.63 

-0.01 

0.37 

case5 

65.85 

3.09 

0.05 

111 

1.11 

4.23 

0.02 

0.36 

case6 

65.85 

5.00 

0.08 

111 

1.18 

6.59 

-0.07 

0.42 

case7 

65.85 

4.87 

0.07 

111 

1.18 

6.61 

0.00 

0.40 

case8 

65.85 

4.66 

0.07 

111 

1.18 

6.42 

0.02 

0.38 

Table  5.10:  Model  skill  parameters  for  the  sensitivity  test  to  surface  elevation  at 
Station  2. 


R 

Drms 

Drms/ R 

V 

Gw 

E'rms 

Lm 

Lrms 

easel 

72.12 

3.11 

0.04 

111 

1.04 

2.62 

0.09 

0.40 

case2 

72.12 

1.89 

0.03 

111 

0.99 

1.84 

-0.05 

0.28 

case3 

72.12 

2.35 

0.03 

111 

1.04 

2.38 

-0.01 

0.28 

case4 

72.12 

2.25 

0.03 

111 

1.02 

1.98 

0.01 

0.28 

case5 

72.12 

2.77 

0.04 

111 

1.03 

2.25 

0.08 

0.39 

case6 

72.12 

3.30 

0.05 

111 

1.03 

2.44 

0.11 

0.42 

case7 

72.12 

3.23 

0.04 

111 

1.04 

2.55 

0.10 

0.41 

case8 

72.12 

3.03 

0.04 

111 

1.04 

2.54 

0.07 

0.38 
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1 : Model  skill  parameters  for  the 

sensi 

tivity 

est  to  u-velocity  at  S 

R 

Drms 

Drms/  R 

V 

Gw 

Drms 

Lm 

Lrms 

easel 

102.2 

12.58 

0.12 

111 

0.84 

10.18 

0.35 

0.59 

case2 

102.2 

12.12 

0.12 

111 

0.87 

8.95 

0.36 

0.60 

case3 

102.2 

7.92 

0.08 

111 

0.97 

5.82 

0.14 

0.42 

case4 

102.2 

7.83 

0.08 

111 
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5.8  Comparison  of  cr-  and  z-baroclinic  approaches 


Recently,  some  problems  on  the  cr  grid  model  have  been  raised.  One  problem 
is  from  the  horizontal  density  gradient  along  the  constant  cr-surface,  which  causes 
a hydrostatic  inconsistency  due  to  the  bottom  topography  effect.  The  hydrostatic 
assumption  in  the  governing  equations  makes  this  problem  more  important.  Another 
problem  is  from  the  advection  terms  in  the  transport  equation,  which  enable  some 
communication  between  the  denser  bottom  water  in  deep  areas  and  the  lighter  bottom 
water  in  adjacent  shallow  areas. 

The  error  caused  by  the  pressure  gradient  force  term  near  steep  topography  in 
the  cr-coordinate  system  commonly  has  been  mentioned.  The  pressure  gradient  force 
term  is  divided  into  two  terms  in  the  cr-coordinate  system,  the  gradient  of  pressure 
along  a constant  a surface  and  the  gradient  of  bottom  topography.  For  example,  the 
horizontal  pressure  gradient  in  the  ^-direction  is  written  in  the  cr-coordinate  system 
as 


dp  dp 

dx  dx 

z 

where  p is  the  pressure,  subscript  z indicates  the  constant  z-coordinate  surface,  H is 
the  total  depth,  a = (z  - f £)///,  and  £ is  the  water  elevation. 

A sufficient  condition  for  a finite  difference  scheme  to  be  hydrostatically  consistent  in 
the  cr-coordinate  system,  is  given  by  Haney  (1990),  as 


cr  dp  dH 


(5.15) 


cr  dH 
~H~dx 


8x  < 8a 


(5.16) 


This  condition,  which  requires  a sufficiently  small  grid  size  for  a given  vertical  grid 
increment,  guarantees  that  the  cr-surface  immediately  below  a given  cr-surface  remains 
below  a given  cr-surface  within  a horizontal  distance  of  one  grid  interval.  If  the  above 


condition  is  not  satisfied,  the  finite  difference  scheme  is  hydrostatically  inconsistent 
and  nonconvergent. 
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To  see  the  problem  associated  with  a grid  model,  I tested  the  effect  of  an  abrupt 
change  of  depth.  Here  a new  terminology  is  defined  : cr-baroclinic  for  the  baroclinic 
term  calculated  along  the  constant  cr-surface  and  z-baroclinic  for  the  baroclinic  term 
calculated  along  the  constant  z-surface.  To  change  cr-baroclinic  to  z-baroclinic,  the 
term  associated  with  depth  was  ignored,  which  is  shown  in  Equation  (5.15).  Now 
density  difference  is  taken  in  the  same  z-level.  Density  at  the  adjacent  point  is 
calculated  through  the  interpolation  of  cr-level.  In  doing  so,  the  horizontal  density 
gradient  is  hydrostatically  consistent. 

Comparison  of  cr-baroclinic  and  z-baroclinic  was  made  with  the  case2  input.  The 
harmonic  constants  of  elevations  at  Station  1,  2,  and  3 show  that  there  is  no  significant 
difference  between  cr-baroclinic  and  z-baroclinic.  The  harmonic  constant  of  velocity 
at  Station  4 also  does  not  show  significant  difference. 

Time  series  of  velocity  and  salinity  fields  over  transect  ‘JB’  were  plotted  in  Fig- 
ures 5.92-5.103  for  three  events  of  cr-baroclinic  with  diffusion,  z-baroclinic  with  dif- 
fusion, and  cr-baroclinic  without  diffusion.  Diffusion  coefficient  used  here  was  10000 
cm2 /sec.  Generally  velocity  in  z-baroclinic  was  larger  than  that  in  cr-baroclinic.  The 
difference  was  large  in  shallow  depth.  The  degree  of  salinity  stratification  was  greater 
in  cr-baroclinic.  This  is  due  to  the  advection  terms  in  transport  equation.  In  the 
cr-grid  system,  this  is  an  evidence  of  communication  between  deep  area  and  shallow 
area.  The  diffusion  effect  was  not  significant  in  the  test  with  diffusion  and  without 
diffusion  in  the  cr-grid  system  . Figure  5.104-5.107  show  salinity  time  series  over  tran- 
sect ‘JB’  at  surface  and  bottom,  comparing  cr-baroclinic  with  diffusion,  z-baroclinic 
with  diffusion,  and  cr-baroclinic  without  diffusion.  To  reduce  the  differences  between 
a-  and  z-surface,  the  condition  of  Equation  (5.16)  should  be  satisfied. 
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Figure  5.92:  Time  series  of  u-velocity  and  salinity  in  z-baroclinic  with  diffusion  at 
grid  point  (65,16). 
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Figure  5.93:  Time  series  of  u-velocity  and  salinity  in  cr-baroclinic  with  diffusion  at 
grid  point  (65,16). 
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Figure  5.94:  Time  series  of  u-velocity  and  salinity  in  cr-baroclinic  without  diffusion 
at  grid  point  (65,16). 
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Figure  5.95:  Time  series  of  u-velocity  and  salinity  in  .z-baroclinic  with  diffusion  at 
grid  point  (65,17). 
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Figure  5.96:  Time  series  of  u-velocity  and  salinity  in  cr-baroclinic  with  diffusion  at 
grid  point  (65,17). 
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Figure  5.97:  Time  series  of  u- velocity  and  salinity  in  cr-baroclinic  without  diffusion 
at  grid  point  (65,17). 
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Figure  5.98:  Time  series  of  u-velocity  and  salinity  in  z-baroclinic  with  diffusion  at 
grid  point  (65,18). 
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Figure  5.99:  Time  series  of  u-velocity  and  salinity  in  cr-baroclinic  with  diffusion  at 
grid  point  (65,18). 
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Figure  5.100:  Time  series  of  u-velocity  and  salinity  in  cr-baroclinic  without  diffusion 
at  grid  point  (65,18). 
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Figure  5.101:  Time  series  of  u-velocity  and  salinity  in  z-baroclinic  with  diffusion  at 
grid  point  (65,19). 
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Figure  5.102:  Time  series  of  u-velocity  and  salinity  in  cr-baroclinic  with  diffusion  at 
grid  point  (65,19). 
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Figure  5.103:  Time  series  of  u-velocity  and  salinity  in  cr-baroclinic  without  diffusion 
at  grid  point  (65,19). 
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Figure  5.104:  Time  series  of  salinity  in  z-baroclinic  with  diffusion,  cr-baroclinic  with 
diffusion,  and  tr-baroclinic  without  diffusion  at  grid  point  (65,16). 
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Figure  5.105:  Time  series  of  salinity  in  z-baroclinic  with  diffusion,  cr-baroclinic  with 
diffusion,  and  cr-baroclinic  without  diffusion  at  grid  point  (65,17). 
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Figure  5.106:  Time  series  of  salinity  in  z-baroclinic  with  diffusion,  cr-baroclinic  with 
diffusion,  and  cr-baroclinic  without  diffusion  at  grid  point  (65,18). 
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Figure  5.107:  Time  series  of  salinity  in  z-baroclinic  with  diffusion,  cr-baroclinic  with 
diffusion,  and  cr-baroclinic  without  diffusion  at  grid  point  (65,19). 


CHAPTER  6 
CONCLUSIONS 


A three-dimensional  curvilinear  hydrodynamic  model  has  been  formulated,  veri- 
fied, and  applied  to  the  James  River,  a partially-mixed  estuary.  The  focus  of  the  cir- 
culation and  salinity  transport  modeling  study  was  on  stratification-destratification 
phenomena  associated  with  the  spring-neap  tidal  cycle.  Additionally,  the  develop- 
ment of  a frontal  system  was  found  near  Newport  News  Point.  The  model  could 
simulate  the  frontal  system. 

The  model  is  based  on  the  time-variable,  three-dimensional  equations  of  conti- 
nuity, momentum,  and  salt  balance.  The  equations  have  been  transformed  using  a 
curvilinear  coordinate  system  in  the  horizontal  and  a cr-stretching  technique  in  the 
vertical.  The  use  of  curvilinear  coordinates  allows  the  user  to  design  the  computa- 
tional mesh  effectively  and  resolve  complex  geometries. 

A model  of  circulation  and  salinity  transport  for  the  partially-mixed  estuary  has 
been  developed,  and  verified  against  analytical  solutions,  and  applied  to  the  James 
River.  Also  the  sensitivity  tests  of  various  parameters  have  been  performed. 

The  model  is  formulated  on  the  basis  of  the  three-dimensional,  time- varying  equa- 
tions of  continuity,  momentum,  and  salt  balance.  The  program  code  can  handle  the 
integration  of  the  governing  equations,  using  variable  weighting  factors  from  the  fully- 
implicit  scheme  to  the  fully-explicit  scheme.  In  doing  so,  the  model  can  be  more 
flexible  to  apply  to  a real  situation.  The  time  step  of  integration  is  more  relaxed  as 
the  weighting  factor  of  implicit  scheme  increases. 

The  model  was  formulated  using  weighting  factor  from  implicit  scheme  (0  = 1) 
to  explicit  scheme  (0  = 0).  With  this  formulation,  the  model  is  more  flexible  when 


252 


253 

applied  to  a special  situation.  In  the  present  study,  the  implicit  integration  scheme 
for  the  momentum  was  used  to  relax  CFL  condition  limiting  the  integration  time 
step.  The  model  solves  the  governing  equations  using  a mode-splitting  technique  in 
which  the  external  or  vertically  averaged  mode  and  the  internal  or  vertical  structure 
are  solved  separately. 

The  vertical  eddy  coefficients  for  momentum  and  salinity  were  determined  by 
the  simplified  second  order  equilibrium  closure  scheme.  So  an  instantaneous  velocity 
was  considered  for  the  eddy  viscosity  and  diffusivity  parameters.  Compared  to  the 
equilibrium  closure  model,  the  kinetic  energy  closure  model  (Sheng  and  Villaret, 
1989)  contains  an  extra  differential  transport  equation  for  q 2 . However,  no  study  has 
unequivocally  demonstrated  that  the  kinetic  energy  closure  model  always  gives  much 
better  results  than  the  equilibrium  closure  model  in  simulating  estuarine  circulation 
and  transport.  And  the  additional  equations  do  not  necessarily  guarantee  better 
model  results. 

A series  of  model  test  runs  was  performed  against  analytic  solutions  to  check  the 
accuracy  of  model  prediction.  Test  of  a forced  constant  surface  slope  showed  that  the 
surface  slope  balanced  the  bottom  stress  in  the  model.  Model  tests  of  seiche,  tidal 
forcing,  wind  forcing,  and  density-induced  forcing  all  showed  good  agreements  with 
the  analytic  solutions. 

Model  calibration  and  validation  was  performed  through  the  quantitative  mea- 
sures of  deviation  between  model  results  and  data.  The  use  of  quantitative  measures  of 
calibration  and  validation  for  hydrodynamic  models  will  enhance  model  performance 
evaluation.  A series  of  sensitivity  tests  were  carried  out  to  determine  the  sensitivity 
to  the  changes  in  some  model  parameters.  Sensitivity  to  bottom  friction  between 
0.002  to  0.004  showed  that  the  differences  were  hardly  noticeable.  With  time-varying 
bottom  friction  calculated  from  the  law  of  the  wall  with  Monin-Oboukhov  length 
scale,  the  model  result  showed  the  best  agreement  with  data.  Roughness  height  was 
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another  important  parameter  to  affect  the  bottom  stress.  With  a small  roughness 
height,  tidal  elevation  increased  and  phase  lag  decreased.  The  bathymetry  was  also 
an  important  factor  to  produce  a good  result,  because  it  is  directly  related  to  the 
bottom  friction.  Space-time  discretization  can  also  result  in  an  improper  simulation 
because  the  discretization  of  the  governing  equations  leads  unavoidably  to  approxi- 
mation errors.  With  small  time  step  and  fine  grid,  the  convergence  and  stability  can 
definitely  be  improved. 

Besides  numerical  constraints,  modeling  of  estuaries  in  three-dimensional  struc- 
ture has  been  hampered  by  lack  of  sufficient  data  for  verification  of  long-term  model 
simulations.  Four  series  of  surveys,  however,  were  available  for  use  of  this  study. 
Each  case  was  conducted  to  make  a vertical-longitudinal  salinity  profile  along  the 
axis  of  the  river.  Application  of  the  model  to  the  spring-neap  stratification  cycle  was 
generally  successful.  The  model  simulation  showed  stratification-destratification  in 
response  to  the  change  of  tidal  range.  The  stratification-destratification  phenomenon 
was  related  to  the  instantaneous  velocity,  which  was  considered  in  the  second  order 
closure  scheme. 

The  longitudinal  motion  of  the  salinity  intrusion  was  simulated  according  to  the 
river  discharge.  The  head  of  salt  intrusion  moved  up  and  down  due  to  the  change 
of  river  discharge.  The  intrusion  moves  downstream  with  inceasing  river  discharge. 
Also  increasing  river  discharge  caused  an  increase  in  relative  stratification. 

Frontal  boundary  was  simulated  in  the  model.  The  front  near  Newport  News 
Point  was  revealed  as  a narrow  zone  of  surface  flow  convergence  and  strong  vertical 
transport. 

It  was  found  that  with  an  estuarine  system  with  a deep  channel  it  is  difficult 
to  predict  the  correct  stratification  degree  and  the  velocity  field  near  the  bottom. 
This  problem  comes  from  grid  resolution  of  channel  depth  and  a transformation. 
Near  steep  topography  the  pressure  gradient  force  is  affected  significantly  with  the 
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gradient  of  bottom  topography  in  cr  coordinate  sytem.  To  reduce  this  problem,  the 
use  of  z-baroclinic  or  z-grid  model  is  highly  recommended. 

Again  hydrodynamical  numerical  modeling  in  partially-mixed  estuaries  is  a very 
complicated  problem.  In  partially  mixed  estuaries,  there  are  freshwater  discharge, 
tidal  forcing,  and  turbulent  kinetic  energy  generated  by  a mean  flow. 
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APPENDIX  A 

HIGHER-ORDER-TERMS  (H.O.T.)  IN  <r-GRID  SYSTEM 


When  second-order  derivative  terms  in  the  momentum  equation  are  transformed 
from  ( x,y,z,t ) to  (x' ,y' ,a,t'),  H.O.T.s  are  generated. 
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APPENDIX  B 

TRANSFORMATION  OF  GOVERNING  EQUATIONS 


B.l  Transformation  Rules 

As  shown  in  Fig.  8 of  Appendix  A of  Sheng  (1986),  the  contravariant  components 
(u‘)  and  physical  components  u(i ) of  the  velocity  vector  in  the  non-Cartesian  system 
are  locally  parallel  or  orthogonal  to  the  grid  lines,  while  the  covariant  components  (u,) 
are  generally  not  parallel  or  orthogonal  to  the  local  grid  lines.  The  three  components 
are  identical  in  a Cartesian  coordinate  system.  The  following  relationships  are  valid 
for  the  three  components  in  a non- Cartesian  system 

u'  = {ga)~l^2u{i)  (no  sum  on  i)  (B.l) 


ui  = (gu)  l/2giMj) 


(no  sum  on  i) 


(B.2) 


u(i)  = Hi 


(B.3) 


where  g is  the  diagonal  element  of  the  metric  tensor  gij 

_ dxl dxn  c 

9ij  = d(i  dz> bfn 

which  for  the  two-dimensional  case  of  interest  is 


(B.4) 
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The  three  components  follow  different  rules  for  transformation  between  the  prototype 
and  the  transformed  plane 


u = 
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dx’u 
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dt’ 

*■  = a?Ui 


(B.7) 


u(0  = (ff*)1/2^7«(i)  (B-8) 

where  the  unbarred  quantities  represent  the  components  in  the  prototype  system  and 
the  barred  quantities  represent  the  components  in  the  transformed  system. 

B.2  Dimensionless  Tensor  Invariant  Governing  Equations 

Before  transforming  the  governing  equations,  it  is  essential  to  first  write  them  in 
tensor-invariant  form,  i.e.,  equations  which  are  independent  of  coordinate  translation 
and  rotation. 

Based  on  extensive  tests,  we  have  concluded  that  using  the  contravariant  com- 
ponents will  yield  the  simplest  form  of  the  governing  equations.  For  simplicity,  we 
will  use  un-barred  quantities  to  denote  the  variables  in  the  transformed  system  unless 
otherwise  indicated. 

Following  the  rules  described  in  the  previous  paragraph,  we  can  obtain  the  fol- 
lowing equations 
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Fvl 

where  d/dxk  is  the  partial  derivative,  gln  is  the  metric  tensor  while  gQ  = J = 
xtDr)  ~ EriUt  is  the  determinant  of  the  metric  tensor,  uk  is  the  contravariant  veloc- 
( ),/  represents  the  covariant  spatial  derivative,  \k  represents  the  contravariant 
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spatial  derivative,  and  e*J  is  the  permutation  tensor  and 
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(B-l  1) 


The  covariant  and  contravariant  differentiations  are  defined  by 
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*'.  = u'  4-  D'  ua 


(B.12) 


S'!1  = gtmS,m 


(B-13) 


where  :j  represents  partial  differentiation  and  D'aj  represents  the  Christoffel  symbol 
of  the  second  kind 

D)k  = ginDnjk  (B.14) 

where  g'n  represents  the  inverse  metric  tensor,  hin,  and  Dnjk  is  the  Christoffel  symbol 
of  the  first  kind 


Dijk  — 2 {.Sij'.k  + 9ik:j  9jk:i ) (B.15) 

B.3  Transformation  of  stress  components 

Using  the  chain  rule  of  partial  differentiation,  the  contravariant  stress  components 
are  transformed  from  the  Cartesian  system  (r1,  r2)  to  the  transformed  system  (r1,  r2), 
and  vice  versa.  The  unbarred  quantities  represent  the  components  in  the  Cartesian 
system  while  the  barred  quantuties  represent  the  components  in  the  transformed 
system. 
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APPENDIX  C 

HORIZONTAL  DIFFUSION  TERMS  IN  BOUNDARY-FITTED  GRIDS 

The  horizontal  diffusion  terms  in  u-velocity  equation  are 
+ 2U12//12  + 1*22^22 

+ (^llu)l^ll  + (DjjU)2i/l2  + {D\2u)\H\2  + (Dj2w)2^22 
T (D\2v)iHu  -f  (D\2v)2Hi2  + (D\2v)iHi2  + ( D\2v)2H22 

— {D\2H\2  + D\2H22)U\ 
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+ \{D\2D22  - D2nD\2)Hn  + (D\xD\2  + D\2D\2  - D\2D\2  - D\2D\2)H12]v 

(C.l) 

The  horizontal  diffusion  terms  in  u-velocity  equation  are 
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where  the  i/tJ  are  inverse  metric  tensor,  the  Dx-k  is  the  Christoffel  symbol,  and  the 
other  subscripts  represent  partial  differentiation. 


APPENDIX  D 

MORE  COMPACT  MOMENTUM  EQUATIONS 


The  u-equation  and  v 


-equation  are  presented  in  the  following 
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+ EfjAu  ■ (Horizontal  Diffusion  of  u ) 
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Kd_  f dv\ 
H2da  { vda] 


+ Eh  Ah  • (Horizontal  Diffusion) 


APPENDIX  E 

FINITE  DIFFERENCE  EQUATIONS  IN  (£,r?,<r)  (THREE-TIME-LEVEL) 


E.l  External  Mode 

The  three-time-level  scheme  of  vertically-integrated  differential  equations  is  writ- 
ten as  follows 
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The  above  equations  are  written  in  finite  difference  forms  as  follows 
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The  above  equations  can  be  written  in  ADI  scheme  as  follows 
X- sweep: 
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f/n+1  is  calculated  again  with  £n+1. 
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Tiv^712^(Cn+1)  + (Vs~oUn+1)  = {yfclT)  + (E.ll) 

where  T\ , T2  and  T3  are  all  constants  between  0 and  1 . For  example,  the  wave  propa- 
gation terms  are  treated  explicitly  if  Tx  = 0,  but  implicitly  if  Tx  = 1.  m can  be  n — 1 
or  n.  Additional  parameters  appearing  in  the  equations  are 
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The  system  of  equations  can  be  rewritten  in  the  following  matrix  form 
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and  F2  and  G2  indicate  all  the  explicit  terms  in  U and  V equations,  respectively, 
and  F3  and  G3  indicate  all  the  explicit  terms  in  the  u and  v equations,  respectively. 
The  subscript  s indicates  a £-point,  subscript  u indicates  a u-point,  and  subscript  v 
indicates  u-point. 
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